1 #include "sico2elmer.h"
2
3 #include "../../config.h"
4
5 /*
6 jv: added fortran function name and char ptr macros to (hopefully) enhance portability
7 */
8
9 /******************************************************
10 Output of SICOPOLIS grid in ELMER POST format
11 *******************************************************/
FC_FUNC(postgrid,POSTGRID)12 void STDCALLBULL FC_FUNC(postgrid,POSTGRID) (float *xi, /* unscaled x coordinate index i: 0,imax */
13 float *eta, /* unscaled y coordinate index j from 0,jmax */
14 float *z_c, /* unscaled z coordinate index i: 0,imax, j: 0,jmax, kc: 0,kcmax */
15 float *z_t, /* unscaled y coordinate index i: 0,imax, j: 0,jmax, kt: 0,kcmax */
16 float *deltaX, /* horizontal grod spacing */
17 int *imax_in, /* grid steps in xi-direction */
18 int *jmax_in, /* grid steps in eta-direction */
19 int *kcmax_in, /* grid steps in z-direction in cold ice layer */
20 int *ktmax_in, /* grid steps in z-direction in temperate ice layer */
21 FC_CHAR_PTR(runname,i1), /*name of run*/
22 FC_CHAR_PTR(ergnum,i2), /*number of file*/
23 int *maske, /*mask of vertex type */
24 int *flag){
25 register int i, j, k, n, element_nr, boundary_nr,kn;
26 int number_of_layers[2], number_of_elements, elements_in_one_layer, number_of_iced_nodes, number_of_nodes[2], nodes_of_element[8], *nodes_of_side_element, number_of_iced_collums, nodes_in_one_layer, number_of_ice_boundaries, *iced, *boundary, *gridmap, auxiliary;
27 int imax, jmax, kcmax, ktmax;
28 float *staggered_grid[2];
29 float actual_scaled_coord[3];
30 char groupid[4], filename[80], yes_no, *suffix=".ep";
31 FILE *ptFile;
32
33 /* constants */
34 imax= *imax_in;
35 jmax= *jmax_in;
36 kcmax= *kcmax_in;
37 ktmax= *ktmax_in;
38 nodes_in_one_layer = (imax+2)*(jmax+2);
39 elements_in_one_layer = (imax+1)*(jmax+1);
40 number_of_layers[0] = kcmax+1;
41 number_of_layers[1] = (*flag)*(ktmax);
42 number_of_nodes[0] = nodes_in_one_layer*number_of_layers[0];
43 number_of_nodes[1] = nodes_in_one_layer*number_of_layers[1];
44 /* print out little summary */
45 printf("---------------------------------------------------------------\n");
46 printf("| Output of SICOPOLIS Grid vor ELMER POST\n");
47 printf("---------------------------------------------------------------\n");
48 printf("| imax/jmax/kcmax/ktmax=%d/%d/%d/%d\n",imax, jmax, kcmax, ktmax);
49 printf("---------------------------------------------------------------\n");
50 printf("| nodes in original grid:\n");
51 printf("| cold layer: %d=%d * %d\n", (imax+1)*(jmax+1)*(kcmax+1), (imax+1)*(jmax+1), (kcmax+1));
52 printf("| temperate layer: %d=%d * %d\n", (imax+1)*(jmax+1)*(ktmax+1), (imax+1)*(jmax+1), (ktmax+1));
53 printf("| -------------\n");
54 printf("| %d=%d * %d\n", (imax+1)*(jmax+1)*(ktmax+2+kcmax), (imax+1)*(jmax+1), ktmax+2+kcmax);
55 printf("---------------------------------------------------------------\n");
56 if (*flag)
57 printf("| Output of temperate layer enabled\n");
58 else printf("| Output of temperate layer disabled\n");
59 printf("---------------------------------------------------------------\n");
60 printf("| nodes in full staggered grid:\n");
61 printf("| cold layer: %d=%d*%d\n", number_of_nodes[0], nodes_in_one_layer, number_of_layers[0]);
62 printf("| temperate layer: %d=%d*%d\n", number_of_nodes[1], nodes_in_one_layer, number_of_layers[1]);
63 printf("| -------------\n");
64 printf("| %d\n", number_of_nodes[1] + number_of_nodes[0]);
65 printf("---------------------------------------------------------------\n");
66 /* allocate memory */
67 staggered_grid[0] = (float *) malloc((size_t) (number_of_nodes[0]*3*sizeof(float)));
68 if (staggered_grid[0] == NULL){
69 printf("ERROR in allocating memory for staggered grid of cold ice layer\n");
70 return;
71 }
72 staggered_grid[1] = (float *) malloc((size_t) (number_of_nodes[1]*3*sizeof(float)));
73 if (staggered_grid[1] == NULL){
74 printf("ERROR in allocating memory for staggered grid of temperate ice layer\n");
75 free(staggered_grid[0]);
76 return;
77 }
78 iced = (int *) malloc((size_t) ((size_t) (imax+1)*(jmax+1)*sizeof(int)));
79 if (iced == NULL){
80 printf("ERROR in allocating memory for glaciation information\n");
81 free(staggered_grid[0]);free(staggered_grid[1]);
82 return;
83 }
84 /* get staggered grid(s) */
85 auxiliary = get_staggered_grid(xi,eta,z_c,imax,jmax,kcmax,deltaX,staggered_grid[0]);
86 /* exit(0); */
87 if (auxiliary!=number_of_nodes[0]){
88 printf("WARNING: number of written %d gridpoints for cold layer does not match calculated %d\n",auxiliary, number_of_nodes[0]);
89 printf("ERROR: Interpolation of staggered grid for cold layer failed!\n");
90 free(staggered_grid[0]);free(staggered_grid[1]);
91 return;
92 }else
93 printf("| succeeded in interpolating staggered grid for cold layer\n");
94 if (*flag){
95 auxiliary = get_staggered_grid(xi,eta,z_t,imax,jmax,(ktmax-1),deltaX,staggered_grid[1]);
96 if(auxiliary!=number_of_nodes[1]){
97 printf("WARNING: number of written %d gridpoints for cold layer does not match calculated %d\n",auxiliary, number_of_nodes[1]);
98 printf("ERROR: Interpolation of staggered grid for temperate layer failed!\n");
99 free(staggered_grid[0]);free(staggered_grid[1]);
100 return;
101 }else
102 printf("| succeeded in interpolating staggered grid for temperate layer\n");
103 }else
104 printf("| no staggered grid for temperate layer interpolated\n");
105 printf("---------------------------------------------------------------\n");
106 /* get glaciation info */
107 number_of_iced_collums=get_glaciation_info(imax,jmax,iced,maske);
108 if (number_of_iced_collums<1){
109 printf("| no glaciation\n");
110 boundary = (int *) malloc((size_t) 4*sizeof(int)); /* dummy array size */
111 /* free(staggered_grid[0]);free(staggered_grid[1]);free(iced); */
112 /* return; */
113 }else{
114 printf("| number of iced colums: %d out of %d (%3.2f %%)\n", number_of_iced_collums, (imax+1)*(jmax+1), ((float) number_of_iced_collums)*100.0/((float) (imax+1)*(jmax+1)));
115 boundary = (int *) malloc((size_t) number_of_iced_collums*4*sizeof(int));
116 }
117 if (boundary==NULL){
118 printf("ERROR in allocation of memory for ice-boundary information\n");
119 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);
120 return;
121 }
122 gridmap = (int *) malloc((size_t) nodes_in_one_layer*sizeof(int));
123 if (gridmap==NULL){
124 printf("ERROR in allocation of memory for staggered grid glaciation information\n");
125 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);
126 return;
127 }
128 number_of_iced_nodes = get_grid_map(imax,jmax,iced,gridmap);
129 if (number_of_iced_nodes>1){
130 /* printf("ERROR in calculation of glaciation information for staggered grid\n"); */
131 /* free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(gridmap); */
132 /* return; */
133 printf("| number of iced nodes in one layer: %d out of %d\n",number_of_iced_nodes, nodes_in_one_layer);
134 }
135
136 number_of_ice_boundaries = get_glaciation_boundary_info(imax,jmax,iced,boundary);
137 if (number_of_ice_boundaries<1){
138 printf("| no glaciation\n");
139 /* printf("ERROR in calculation of boundaries\n"); */
140 /* free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap); */
141 /* return; */
142 }else{
143 printf("| number of boundary elements lines: %d (%f km)\n", number_of_ice_boundaries, ((float) number_of_ice_boundaries)*(*deltaX));
144 }
145 /* calculate constants for output mesh*/
146 if (*flag){
147 number_of_nodes[0] = number_of_iced_nodes*(kcmax+1);
148 number_of_nodes[1] = number_of_iced_nodes*(ktmax-1)+nodes_in_one_layer;
149 }else{
150 number_of_nodes[0] = number_of_iced_nodes*kcmax;
151 number_of_nodes[1] = nodes_in_one_layer;
152 }
153 number_of_elements = elements_in_one_layer + (1+(*flag))*number_of_iced_collums + (kcmax+(*flag)*ktmax)*number_of_iced_collums + number_of_ice_boundaries*(kcmax+(*flag)*ktmax);
154 nodes_of_side_element = (int *) malloc((size_t) number_of_ice_boundaries*4*(kcmax+(*flag)*ktmax)*sizeof(int));
155 if (nodes_of_side_element == NULL){
156 printf("ERROR in allocation of side element array\n");
157 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);
158 return;
159 }
160 printf("---------------------------------------------------------------\n");
161 printf("| number of nodes:\n");
162 printf("| in cold layer %d\n",number_of_nodes[0]+(1-(*flag))*nodes_in_one_layer);
163 printf("| in temperate layer %d\n",number_of_nodes[1]+(*flag)*nodes_in_one_layer);
164 printf("| -----------\n");
165 printf("| %d\n",number_of_nodes[0]+number_of_nodes[1]);
166 printf("---------------------------------------------------------------\n");
167 printf("| number of elements:\n");
168 printf("| in cold layer volume %d\n",kcmax*number_of_iced_collums);
169 printf("| in temperate layer volume %d\n",(*flag)*ktmax*number_of_iced_collums);
170 printf("| at base %d\n",elements_in_one_layer);
171 printf("| at CTS %d\n",(*flag)*number_of_iced_collums);
172 printf("| at free surface %d\n",number_of_iced_collums);
173 printf("| on boundary of cold layer %d\n",number_of_ice_boundaries*(kcmax));
174 printf("| on boundary of temperate layer %d\n",(*flag)*number_of_ice_boundaries*(ktmax));
175 printf("| -----------\n");
176 printf("| %d\n",number_of_elements);
177 printf("---------------------------------------------------------------\n");
178 /* write mesh file header*/
179 sprintf(filename,"%5s%2s%s", runname, ergnum, suffix);
180 printf("| Writing mesh-file %s.\n",filename);
181 if((ptFile=fopen(filename, "w"))==NULL){
182 printf("\a Could not open Elmer mesh-file file for writing!\n");
183 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);
184 return;
185 }
186 fprintf(ptFile, "%d %d 1 1\n", number_of_nodes[0]+number_of_nodes[1], number_of_elements);
187 /* write nodes */
188 if (*flag){
189 for(j=0,n=0;j<jmax+2;++j){
190 for (i=0;i<imax+2;++i){
191 fprintf(ptFile, "%6.4e %6.4e %6.4e\n", staggered_grid[1][(j*(imax+2) + i)*3 + 0], staggered_grid[1][(j*(imax+2) + i)*3 + 1], staggered_grid[1][(j*(imax+2) + i)*3 + 2]);++n;
192 }
193 }
194 for (k=1;k<ktmax;++k){
195 for(j=0;j<jmax+2;++j){
196 for (i=0;i<imax+2;++i){
197 if (gridmap[j*(imax+2) + i]!=-1){
198 fprintf(ptFile, "%6.4e %6.4e %6.4e\n", staggered_grid[1][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 0], staggered_grid[1][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 1], staggered_grid[1][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 2]);
199 ++n;
200 }
201 }
202 }
203 }
204 for (k=0;k<kcmax+1;++k){
205 for(j=0;j<jmax+2;++j){
206 for (i=0;i<imax+2;++i){
207 if (gridmap[j*(imax+2) + i]!=-1){
208 fprintf(ptFile, "%6.4e %6.4e %6.4e\n", staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 0], staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 1], staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 2]);
209 ++n;
210 }
211 }
212 }
213 }
214 }else{
215 for(j=0,n=0;j<jmax+2;++j){
216 for (i=0;i<imax+2;++i){
217 fprintf(ptFile, "%6.4e %6.4e %6.4e\n", staggered_grid[0][(j*(imax+2) + i)*3 + 0], staggered_grid[0][(j*(imax+2) + i)*3 + 1], staggered_grid[0][(j*(imax+2) + i)*3 + 2]);++n;
218 }
219 }
220 for (k=1;k<kcmax+1;++k){
221 for(j=0;j<jmax+2;++j){
222 for (i=0;i<imax+2;++i){
223 if (gridmap[j*(imax+2) + i]!=-1){
224 fprintf(ptFile, "%6.4e %6.4e %6.4e\n", staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 0], staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 1], staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 2]);
225 ++n;
226 }
227 }
228 }
229 }
230 }
231 if (n!=number_of_nodes[0]+number_of_nodes[1]){
232 printf("WARNING: Number of written nodes %d does not match calculated %d\n", n, number_of_nodes[0]+number_of_nodes[1]);
233 fclose(ptFile);
234 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);
235 return;
236 }else printf("| %d nodes written\n",n);
237 /* write elements */
238 /******** Order of Elements *****
239 * 7 6 *
240 * +---------+ E=0 *
241 * |\ \ N=1 *
242 * | \ | \ W=2 *
243 * j | \4 \5 S=3 *
244 * ^ | +-----+---+ *
245 * \ | | | *
246 * \| |N | | *
247 * + - + - -S+ | *
248 * 3\ | 2\ | *
249 * \ | E| *
250 * W \| \| *
251 * +---------+ ����>i *
252 * 0 S 1 *
253 ********************************/
254 boundary_nr=0;
255 element_nr=0;
256 if (*flag){
257 sprintf(groupid,"temp");
258 for (j=0;j<jmax+1;++j){
259 for (i=0;i<imax+1;++i){
260 if (iced[j*(imax+1)+i]!=-1){
261 nodes_of_element[0] = j*(imax+2)+i;
262 nodes_of_element[1] = j*(imax+2)+i+1;
263 nodes_of_element[2] = (j+1)*(imax+2)+i+1;
264 nodes_of_element[3] = (j+1)*(imax+2)+i;
265 nodes_of_element[4] = nodes_in_one_layer + gridmap[j*(imax+2)+i];
266 nodes_of_element[5] = nodes_in_one_layer + gridmap[j*(imax+2)+i+1];
267 nodes_of_element[6] = nodes_in_one_layer + gridmap[(j+1)*(imax+2)+i+1];
268 nodes_of_element[7] = nodes_in_one_layer + gridmap[(j+1)*(imax+2)+i];
269 fprintf(ptFile,"%s 808 %d %d %d %d %d %d %d %d\n", groupid,
270 nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3],
271 nodes_of_element[4], nodes_of_element[5], nodes_of_element[6], nodes_of_element[7]);
272 ++element_nr;
273 if (boundary[iced[(j*(imax+1)+i)]*4+0]){/* east */
274 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[2];
275 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[1];
276 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[5];
277 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[6];
278 ++boundary_nr;
279 }
280 if (boundary[iced[(j*(imax+1)+i)]*4+1]){/* north */
281 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
282 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[2];
283 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[6];
284 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
285 ++boundary_nr;
286 }
287 if (boundary[iced[(j*(imax+1)+i)]*4+2]){/* west */
288 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
289 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
290 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
291 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
292 ++boundary_nr;
293 }
294 if (boundary[iced[(j*(imax+1)+i)]*4+3]){/* south */
295 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[1];
296 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
297 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
298 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[5];
299 ++boundary_nr;
300 }
301 }
302 }
303 }
304 for (k=1;k<ktmax;++k){
305 for (j=0;j<jmax+1;++j){
306 for (i=0;i<imax+1;++i){
307 if (iced[j*(imax+1)+i]!=-1){
308 for (n=0; n<2; n++){/* lower level: n=0; upper level n=1; each counterclkws */
309 nodes_of_element[n*4] = nodes_in_one_layer + (k-1+n)*number_of_iced_nodes + gridmap[j*(imax+2)+i];
310 nodes_of_element[n*4+1] = nodes_in_one_layer + (k-1+n)*number_of_iced_nodes + gridmap[j*(imax+2)+i+1];
311 nodes_of_element[n*4+2] = nodes_in_one_layer + (k-1+n)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i+1];
312 nodes_of_element[n*4+3] = nodes_in_one_layer + (k-1+n)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i];
313 }
314 fprintf(ptFile,"%s 808 %d %d %d %d %d %d %d %d\n", groupid,
315 nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3],
316 nodes_of_element[4], nodes_of_element[5], nodes_of_element[6], nodes_of_element[7]);
317 ++element_nr;
318 if (boundary[iced[(j*(imax+1)+i)]*4+0]){/* east */
319 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[2];
320 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[1];
321 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[5];
322 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[6];
323 ++boundary_nr;
324 }
325 if (boundary[iced[(j*(imax+1)+i)]*4+1]){/* north */
326 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
327 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[2];
328 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[6];
329 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
330 ++boundary_nr;
331 }
332 if (boundary[iced[(j*(imax+1)+i)]*4+2]){/* west */
333 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
334 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
335 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
336 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
337 ++boundary_nr;
338 }
339 if (boundary[iced[(j*(imax+1)+i)]*4+3]){/* south */
340 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[1];
341 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
342 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
343 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[5];
344 ++boundary_nr;
345 }
346 }
347 }
348 }
349 }
350 }
351 sprintf(groupid,"cold");
352 if (*flag) kn=0;
353 else{
354 for (j=0;j<jmax+1;++j){
355 for (i=0;i<imax+1;++i){
356 if (iced[j*(imax+1)+i]!=-1){
357 nodes_of_element[0] = j*(imax+2)+i;
358 nodes_of_element[1] = j*(imax+2)+i+1;
359 nodes_of_element[2] = (j+1)*(imax+2)+i+1;
360 nodes_of_element[3] = (j+1)*(imax+2)+i;
361 nodes_of_element[4] = nodes_in_one_layer + gridmap[j*(imax+2)+i];
362 nodes_of_element[5] = nodes_in_one_layer + gridmap[j*(imax+2)+i+1];
363 nodes_of_element[6] = nodes_in_one_layer + gridmap[(j+1)*(imax+2)+i+1];
364 nodes_of_element[7] = nodes_in_one_layer + gridmap[(j+1)*(imax+2)+i];
365 fprintf(ptFile,"%s 808 %d %d %d %d %d %d %d %d\n", groupid,
366 nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3],
367 nodes_of_element[4], nodes_of_element[5], nodes_of_element[6], nodes_of_element[7]);
368 ++element_nr;
369 if (boundary[iced[(j*(imax+1)+i)]*4+0]){/* east */
370 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[2];
371 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[1];
372 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[5];
373 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[6];
374 ++boundary_nr;
375 }
376 if (boundary[iced[(j*(imax+1)+i)]*4+1]){/* north */
377 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
378 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[2];
379 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[6];
380 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
381 ++boundary_nr;
382 }
383 if (boundary[iced[(j*(imax+1)+i)]*4+2]){/* west */
384 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
385 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
386 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
387 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
388 ++boundary_nr;
389 }
390 if (boundary[iced[(j*(imax+1)+i)]*4+3]){/* south */
391 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[1];
392 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
393 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
394 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[5];
395 ++boundary_nr;
396 }
397 }
398 }
399 }
400 kn=1;
401 }
402 for (k=kn;k<kcmax;++k){
403 for (j=0;j<jmax+1;++j){
404 for (i=0;i<imax+1;++i){
405 if (iced[j*(imax+1)+i]!=-1){
406 for (n=0; n<2; n++){/* lower level: n=0; upper level n=1; each counterclkws */
407 nodes_of_element[n*4] = number_of_nodes[1] + (k-kn+n)*number_of_iced_nodes + gridmap[j*(imax+2)+i];
408 nodes_of_element[n*4+1] = number_of_nodes[1] + (k-kn+n)*number_of_iced_nodes + gridmap[j*(imax+2)+i+1];
409 nodes_of_element[n*4+2] = number_of_nodes[1] + (k-kn+n)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i+1];
410 nodes_of_element[n*4+3] = number_of_nodes[1] + (k-kn+n)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i];
411 }
412 fprintf(ptFile,"%s 808 %d %d %d %d %d %d %d %d\n", groupid,
413 nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3],
414 nodes_of_element[4], nodes_of_element[5], nodes_of_element[6], nodes_of_element[7]);
415 ++element_nr;
416 if (boundary[iced[(j*(imax+1)+i)]*4+0]){/* east */
417 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[2];
418 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[1];
419 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[5];
420 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[6];
421 ++boundary_nr;
422 }
423 if (boundary[iced[(j*(imax+1)+i)]*4+1]){/* north */
424 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
425 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[2];
426 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[6];
427 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
428 ++boundary_nr;
429 }
430 if (boundary[iced[(j*(imax+1)+i)]*4+2]){/* west */
431 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
432 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
433 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
434 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
435 ++boundary_nr;
436 }
437 if (boundary[iced[(j*(imax+1)+i)]*4+3]){/* south */
438 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[1];
439 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
440 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
441 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[5];
442 ++boundary_nr;
443 }
444 }
445 }
446 }
447 }
448 sprintf(groupid,"base");
449 for (j=0;j<jmax+1;++j){
450 for (i=0;i<imax+1;++i){
451 nodes_of_element[0]=j*(imax+2)+i;
452 nodes_of_element[1]=j*(imax+2)+i+1;
453 nodes_of_element[2]=(j+1)*(imax+2)+i+1;
454 nodes_of_element[3]=(j+1)*(imax+2)+i;
455 fprintf(ptFile,"%s 404 %d %d %d %d\n", groupid, nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3]);
456 ++element_nr;
457 }
458 }
459 if (*flag){
460 sprintf(groupid,"cts");
461 for (j=0;j<jmax+1;++j){
462 for (i=0;i<imax+1;++i){
463 if (iced[j*(imax+1)+i]!=-1){
464 nodes_of_element[0] = number_of_nodes[1] + gridmap[j*(imax+2)+i];
465 nodes_of_element[1] = number_of_nodes[1] + gridmap[j*(imax+2)+i+1];
466 nodes_of_element[2] = number_of_nodes[1] + gridmap[(j+1)*(imax+2)+i+1];
467 nodes_of_element[3] = number_of_nodes[1] + gridmap[(j+1)*(imax+2)+i];
468 fprintf(ptFile,"%s 404 %d %d %d %d\n", groupid, nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3]);
469 ++element_nr;
470 }
471 }
472 }
473 }
474 sprintf(groupid,"free");
475 for (j=0;j<jmax+1;++j){
476 for (i=0;i<imax+1;++i){
477 if (iced[j*(imax+1)+i]!=-1){
478 nodes_of_element[0] = number_of_nodes[1] + (kcmax)*number_of_iced_nodes + gridmap[j*(imax+2)+i];
479 nodes_of_element[1] = number_of_nodes[1] + (kcmax)*number_of_iced_nodes + gridmap[j*(imax+2)+i+1];
480 nodes_of_element[2] = number_of_nodes[1] + (kcmax)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i+1];
481 nodes_of_element[3] = number_of_nodes[1] + (kcmax)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i];
482 fprintf(ptFile,"%s 404 %d %d %d %d\n", groupid, nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3]);
483 ++element_nr;
484 }
485 }
486 }
487 if (*flag){
488 sprintf(groupid,"side_t");
489 for (n=0;n<number_of_ice_boundaries*ktmax;++n){
490 fprintf(ptFile,"%s 404 %d %d %d %d\n", groupid, nodes_of_side_element[n*4+0], nodes_of_side_element[n*4+1], nodes_of_side_element[n*4+2], nodes_of_side_element[n*4+3]);
491 ++element_nr;
492 }
493 sprintf(groupid,"side_c");
494 for (n=number_of_ice_boundaries*ktmax;n<number_of_ice_boundaries*(kcmax+ktmax);++n){
495 fprintf(ptFile,"%s 404 %d %d %d %d\n", groupid, nodes_of_side_element[n*4+0], nodes_of_side_element[n*4+1], nodes_of_side_element[n*4+2], nodes_of_side_element[n*4+3]);
496 ++element_nr;
497 }
498 }else{
499 sprintf(groupid,"side_c");
500 for (n=0;n<number_of_ice_boundaries*kcmax;++n){
501 fprintf(ptFile,"%s 404 %d %d %d %d\n", groupid, nodes_of_side_element[n*4+0], nodes_of_side_element[n*4+1], nodes_of_side_element[n*4+2], nodes_of_side_element[n*4+3]);
502 ++element_nr;
503 }
504 }
505 if (element_nr!=number_of_elements)
506 printf("WARNING: number of written elements %d does not match calculated %d\n", element_nr, number_of_elements);
507 else printf("| %d elements written\n", element_nr);
508 printf("---------------------------------------------------------------\n");
509 fclose(ptFile);
510 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);
511 return;
512 }
513 /*******************************************************
514 write data for ELMER post
515 ********************************************************/
FC_FUNC(elmerdata,ELMERDATA)516 void STDCALLBULL FC_FUNC(elmerdata,ELMERDATA) (int *imax_in, /* grid steps in xi-direction */
517 int *jmax_in, /* grid steps in eta-direction */
518 int *kcmax_in, /* grid steps in z-direction in cold ice layer */
519 int *ktmax_in, /* grid steps in z-direction in temperate ice layer */
520 float *z_c, /* z-coordinate in cold layer for given i,j-position in plane and kc in vertical*/
521 float *z_t, /* z-coordinate in temperate layer for given i,j-position in plane and kc in vertical*/
522 float *vx_c, /* velocity component in x for cold region for given i,j-position in plane and kc in vertical */
523 float *vy_c, /* velocity component in y for cold region for given i,j-position in plane and kc in vertical */
524 float *vz_c, /* velocity component in z for cold region for given i,j-position in plane and kc in vertical */
525 float *age_c, /* age in cold region for given i,j-position in plane and kc in vertical */
526 float *temp_c, /* temperature in cold region for given i,j-position in plane and kt in vertical */
527 float *vx_t, /* velocity component in x for temperated region for given i,j-position in plane and kt in vertical */
528 float *vy_t, /* velocity component in y for temperated region for given i,j-position in plane and kt in vertical */
529 float *vz_t, /* velocity component in z for temperated region for given i,j-position in plane and kt in vertical */
530 float *temp_t_m, /* melting temperature in temperate ice region for given i,j-position in plane and kt in vertical */
531 float *age_t, /* age in temperate region for given i,j-position in plane and kc in vertical */
532 float *omega_t, /* H2O mass fraction in temperated region for given i,j-position in plane and kt in vertical */
533 float *Q_bm, /* production rate of melting-water at bottom for given i,j-position in plane */
534 float *Q_tld, /* water drainage rate from the temperated region for given i,j-position in plane */
535 float *am_perp, /* ice volume flux through CTS for given i,j-position in plane */
536 float *qx, /* mass-flux in x-direction for given i,j-position in plane */
537 float *qy, /* mass-flux in y-direction for given i,j-position in plane */
538 int *n_cts, /* polythermal conditions for given i,j-position at base (-1 cold ice; 0 temp. ice base with cold ice above; 1 temp. ice base with temperate ice above; */
539 int *maske, /* glaciation information for given i,j-position at base (glaciated=0, ice-free=1, 2=sea-point, 3=floating ice) */
540 FC_CHAR_PTR(runname,runname_l), /*name of run*/
541 FC_CHAR_PTR(ergnum,ergnum_l), /*number of file*/
542 int *flag){
543 register int i, j, k, n;
544 int imax, jmax, kcmax, ktmax, kcmin, offset, nodes_in_temp, nodes_in_cold, elements_in_layer, number_of_iced_nodes, number_of_written_nodes, *gridmap, ok;
545 int nodes_in_one_layer, number_of_nodes, nodes_in_layer_of_staggered_grid, number_of_iced_collums, *iced, auxiliary, number_of_properties;
546 float *array_for_output, *array_for_interpolation, *age, *temperature, *flux[2], *velocity[3], *height, *drainage, *melt, *ice_land_sea, *type_of_base, *float_cts, *float_maske;
547 char data_filename[80], *suffix=".dat";
548 FILE *ptFile;
549 /* constants */
550 imax=*imax_in;
551 jmax=*jmax_in;
552 kcmax=*kcmax_in;
553 ktmax=*ktmax_in;
554 elements_in_layer=(imax+1)*(jmax+1);
555 nodes_in_layer_of_staggered_grid = (imax+2)*(jmax+2);
556 nodes_in_temp = (*flag)*nodes_in_layer_of_staggered_grid*ktmax;
557 nodes_in_cold = nodes_in_layer_of_staggered_grid*(kcmax+1);
558 number_of_nodes = nodes_in_temp + nodes_in_cold;
559 number_of_properties= (NVEC2*2 + NVEC3*3 + NSCAL);
560 /* print out little summary */
561 printf("---------------------------------------------------------------\n");
562 printf("| Output of SICOPOLIS Data for ELMER POST\n");
563 printf("---------------------------------------------------------------\n");
564 printf("| imax/jmax/kcmax/ktmax=%d/%d/%d/%d\n",imax, jmax, kcmax, ktmax);
565 printf("---------------------------------------------------------------\n");
566 printf("| nodes in original grid:\n");
567 printf("| cold layer: %d=%d * %d\n", (imax+1)*(jmax+1)*(kcmax+1), (imax+1)*(jmax+1), (kcmax+1));
568 printf("| temperate layer: %d=%d * %d\n", (imax+1)*(jmax+1)*(ktmax+1), (imax+1)*(jmax+1), (ktmax+1));
569 printf("| -------------\n");
570 printf("| %d=%d * %d\n", (imax+1)*(jmax+1)*(ktmax+2+kcmax), (imax+1)*(jmax+1), ktmax+2+kcmax);
571 printf("---------------------------------------------------------------\n");
572 if (*flag)
573 printf("| Output of temperate layer enabled\n");
574 else printf("| Output of temperate layer disabled\n");
575 printf("---------------------------------------------------------------\n");
576 printf("| nodes in full staggered grid:\n");
577 printf("| cold layer: %d\n", nodes_in_cold);
578 printf("| temperate layer: %d\n", nodes_in_temp);
579 printf("| -------------\n");
580 printf("| %d\n",number_of_nodes);
581 printf("---------------------------------------------------------------\n");
582 if ((float_maske = (float *) malloc((size_t) nodes_in_layer_of_staggered_grid*sizeof(float)))==NULL){
583 printf("ERROR in allocation of memory!\n");
584 return;
585 }
586 if((float_cts = (float *) malloc((size_t) nodes_in_layer_of_staggered_grid*sizeof(float)))==NULL){
587 printf("ERROR in allocation of memory!\n");
588 free(float_maske);
589 return;
590 }
591 if ((array_for_interpolation = (float *) malloc((size_t) number_of_nodes*number_of_properties*sizeof(float)))==NULL){
592 printf("ERROR in allocation of memory for interpolating data on staggered grid!\n");
593 free(float_cts);
594 free(float_maske);
595 return;
596 }
597 if ((iced = (int *) malloc((size_t) nodes_in_layer_of_staggered_grid*sizeof(int)))==NULL){
598 printf("ERROR in allocation of memory!\n");
599 free(float_cts);
600 free(float_maske);
601 free(array_for_interpolation);
602 }
603 gridmap = (int *) malloc((size_t) nodes_in_layer_of_staggered_grid*sizeof(int));
604 if (gridmap==NULL){
605 printf("ERROR in allocation of memory!\n");
606 free(float_cts);
607 free(float_maske);
608 free(array_for_interpolation);
609 free(iced);
610 }
611 number_of_iced_collums=get_glaciation_info(imax,jmax,iced,maske);
612 number_of_iced_nodes = get_grid_map(imax,jmax,iced,gridmap);
613 if (*flag) number_of_written_nodes=nodes_in_layer_of_staggered_grid+number_of_iced_nodes*(ktmax+kcmax);
614 else number_of_written_nodes=nodes_in_layer_of_staggered_grid+number_of_iced_nodes*(kcmax);
615 printf("| number of iced colums: %d out of %d (%3.2f %%)\n", number_of_iced_collums, (imax+1)*(jmax+1), ((float) number_of_iced_collums)*100.0/((float) (imax+1)*(jmax+1)));
616 printf("| number of iced nodes in layer of staggered grid: %d out of %d\n", number_of_iced_nodes, nodes_in_layer_of_staggered_grid);
617 printf("| number of nodes in written grid: %d of %d\n", number_of_written_nodes, number_of_nodes);
618 printf("| number of properties written for each node: %d\n", number_of_properties);
619 printf("| size of output array =%d=%d*%d\n", number_of_written_nodes*number_of_properties, number_of_written_nodes,number_of_properties);
620 printf("---------------------------------------------------------------\n");
621 array_for_output = (float *) malloc((size_t) number_of_written_nodes*number_of_properties*sizeof(float));
622 if (array_for_output==NULL){
623 printf("ERROR in allocation of memory for output of data!\n");
624 free(float_cts);
625 free(float_maske);
626 free(array_for_interpolation);
627 free(iced);
628 free(gridmap);
629 return;
630 }
631 /* assign pointers to array of output */
632 height = (float *) &array_for_interpolation[0];
633 velocity[0] = (float *) &array_for_interpolation[number_of_nodes];
634 velocity[1] = (float *) &array_for_interpolation[2*number_of_nodes];
635 velocity[2] = (float *) &array_for_interpolation[3*number_of_nodes];
636 drainage = (float *) &array_for_interpolation[4*number_of_nodes];
637 melt = (float *) &array_for_interpolation[5*number_of_nodes];
638 ice_land_sea = (float *) &array_for_interpolation[6*number_of_nodes];
639 type_of_base = (float *) &array_for_interpolation[7*number_of_nodes];
640 temperature = (float *) &array_for_interpolation[8*number_of_nodes];
641 age = (float *) &array_for_interpolation[9*number_of_nodes];
642 flux[0] = (float *) &array_for_interpolation[10*number_of_nodes];
643 flux[1] = (float *) &array_for_interpolation[11*number_of_nodes];
644
645 make_float_from_integer_scalar_field(maske, float_maske, nodes_in_layer_of_staggered_grid, 1);
646 make_float_from_integer_scalar_field(n_cts, float_cts, nodes_in_layer_of_staggered_grid, 0);
647
648 /* set not used 3d parts of 2d arrays to default value */
649 for (n=nodes_in_layer_of_staggered_grid; n<number_of_nodes;++n){
650 ice_land_sea[n]=0.0;
651 type_of_base[n]=2.0;
652 drainage[n]=0.0;
653 melt[n] = 0.0;
654 flux[0][n]=0.0;
655 flux[1][n]=0.0;
656 }
657 /* interpolate properties on staggered grid */
658 ok=1;
659 if (*flag){
660 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, ktmax-1, z_t, height);
661 if (auxiliary!=nodes_in_temp){
662 printf("WARNING: number of returned values %d of interpolated height in temperate layer does not match calculated %d\n",auxiliary, nodes_in_temp);
663 ok=0;
664 }
665 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, ktmax-1, vx_t, &velocity[0][0]);
666 if (auxiliary!=nodes_in_temp){
667 printf("WARNING: number of returned values %d of interpolated velocity_x in temperate layer does not match calculated %d\n", auxiliary, nodes_in_temp);
668 ok=0;
669 }
670 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, ktmax-1, vy_t, &velocity[1][0]);
671 if (auxiliary!=nodes_in_temp){
672 printf("WARNING: number of returned values %d of interpolated velocity_y in temperate layer does not match calculated %d\n",auxiliary, nodes_in_temp);
673 ok=0;
674 }
675 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, ktmax-1, vz_t, &velocity[2][0]);
676 if (auxiliary!=nodes_in_temp){
677 printf("WARNING: number of returned values %d of interpolated velocity_z in temperate layer does not match calculated %d\n",auxiliary, nodes_in_temp);
678 ok=0;
679 }
680 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, ktmax-1, temp_t_m, temperature);
681 if (auxiliary!=nodes_in_temp){
682 printf("WARNING: number of returned values %d of interpolated temperature in temperate layer does not match calculated %d\n",auxiliary, nodes_in_temp);
683 ok=0;
684 }
685 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, ktmax-1, age_t, age);
686 if (auxiliary!=nodes_in_temp){
687 printf("WARNING: number of returned values %d of interpolated age of ice in temperate layer does not match calculated %d\n",auxiliary, nodes_in_temp);
688 ok=0;
689 }
690 }
691 if (!(ok)){
692 printf("ERROR in interpolating data on staggered grid in temperate layer!\n");
693 free(array_for_interpolation);
694 free(float_maske);
695 free(float_cts);
696 free(iced);
697 free(gridmap);
698 free(array_for_output);
699 return;
700 }
701 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, 0, Q_tld, drainage);
702 if (auxiliary!=nodes_in_layer_of_staggered_grid){
703 printf("WARNING: number of returned values %d of interpolated drainage at base layer does not match calculated %d\n",auxiliary, nodes_in_layer_of_staggered_grid);
704 ok=0;
705 }
706 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, 0, Q_bm, melt);
707 if (auxiliary!=nodes_in_layer_of_staggered_grid){
708 printf("WARNING: number of returned values %d of interpolated melting rate at base layer does not match calculated %d\n",auxiliary, nodes_in_layer_of_staggered_grid);
709 ok=0;
710 }
711 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, 0, float_maske, ice_land_sea);
712 if (auxiliary!=nodes_in_layer_of_staggered_grid){
713 printf("WARNING: number of returned values %d of interpolated ice-sea-land mask at base does not match calculated %d\n",auxiliary, nodes_in_layer_of_staggered_grid);
714 ok=0;
715 }
716 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, 0, float_cts, type_of_base);
717 if (auxiliary!=nodes_in_layer_of_staggered_grid){
718 printf("WARNING: number of returned values %d of interpolated ype of base does not match calculated %d\n",auxiliary, nodes_in_layer_of_staggered_grid);
719 ok=0;
720 }
721 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, 0, qx, flux[0]);
722 if (auxiliary!=nodes_in_layer_of_staggered_grid){
723 printf("WARNING: number of returned values %d of interpolated flux_x at base does not match calculated %d\n",auxiliary, nodes_in_layer_of_staggered_grid);
724 ok=0;
725 }
726 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, 0, qy, flux[1]);
727 if (auxiliary!=nodes_in_layer_of_staggered_grid){
728 printf("WARNING: number of returned values %d of interpolated flux_y at base does not match calculated %d\n",auxiliary, nodes_in_layer_of_staggered_grid);
729 ok=0;
730 }
731 if (!(ok)){
732 printf("ERROR in interpolating data on staggered grid at base!\n");
733 free(array_for_interpolation);
734 free(float_maske);
735 free(float_cts);
736 free(iced);
737 free(gridmap);
738 free(array_for_output);
739 return;
740 }
741 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, kcmax, z_c, &height[nodes_in_temp]);
742 if (auxiliary!=nodes_in_cold){
743 printf("WARNING: number of returned values %d of interpolated height in cold layer does not match calculated %d\n",auxiliary, nodes_in_cold);
744 ok=0;
745 }
746 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, kcmax, vx_c, &velocity[0][nodes_in_temp]);
747 if (auxiliary!=nodes_in_cold){
748 printf("WARNING: number of returned values %d of interpolated velocity_x in cold layer does not match calculated %d\n",auxiliary, nodes_in_cold);
749 ok=0;
750 }
751 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, kcmax, vy_c, &velocity[1][nodes_in_temp]);
752 if (auxiliary!=nodes_in_cold){
753 printf("WARNING: number of returned values %d of interpolated velocity_y in cold layer does not match calculated %d\n",auxiliary, nodes_in_cold);
754 ok=0;
755 }
756 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, kcmax, vz_c, &velocity[2][nodes_in_temp]);
757 if (auxiliary!=nodes_in_cold){
758 printf("WARNING: number of returned values %d of interpolated velocity_z in cold layer does not match calculated %d\n",auxiliary, nodes_in_cold);
759 ok=0;
760 }
761
762 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, kcmax, temp_c, &temperature[nodes_in_temp]);
763 if (auxiliary!=nodes_in_cold){
764 printf("WARNING: number of returned values %d of interpolated temperature in cold layer does not match calculated %d\n",auxiliary, nodes_in_cold);
765 ok=0;
766 }
767 auxiliary=get_interpolated_property_on_staggered_grid(imax, jmax, kcmax, age_c, &age[nodes_in_temp]);
768 if (auxiliary!=nodes_in_cold){
769 printf("WARNING: number of returned values %d of age of ice interpolated in cold layer does not match calculated %d\n",auxiliary, nodes_in_cold);
770 ok=0;
771 }
772 if (!(ok)){
773 printf("ERROR in interpolating data on staggered grid in cold layer!\n");
774 free(array_for_interpolation);
775 free(float_maske);
776 free(float_cts);
777 free(iced);
778 free(gridmap);
779 free(array_for_output);
780 return;
781 }
782 /* write to output array */
783 for(j=0,n=0;j<jmax+2;++j){
784 for(i=0;i<imax+2;++i){
785 array_for_output[(j*(imax+2)+i)*number_of_properties+0]=height[j*(imax+2)+i];
786 array_for_output[(j*(imax+2)+i)*number_of_properties+1]=velocity[0][j*(imax+2)+i];
787 array_for_output[(j*(imax+2)+i)*number_of_properties+2]=velocity[1][j*(imax+2)+i];
788 array_for_output[(j*(imax+2)+i)*number_of_properties+3]=velocity[2][j*(imax+2)+i];
789 array_for_output[(j*(imax+2)+i)*number_of_properties+4]=drainage[j*(imax+2)+i];
790 array_for_output[(j*(imax+2)+i)*number_of_properties+5]=melt[j*(imax+2)+i];
791 array_for_output[(j*(imax+2)+i)*number_of_properties+6]=ice_land_sea[j*(imax+2)+i];
792 array_for_output[(j*(imax+2)+i)*number_of_properties+7]=type_of_base[j*(imax+2)+i];
793 array_for_output[(j*(imax+2)+i)*number_of_properties+8]=temperature[j*(imax+2)+i];
794 array_for_output[(j*(imax+2)+i)*number_of_properties+9]=age[j*(imax+2)+i];
795 array_for_output[(j*(imax+2)+i)*number_of_properties+10]=flux[0][j*(imax+2)+i];
796 array_for_output[(j*(imax+2)+i)*number_of_properties+11]=flux[1][j*(imax+2)+i];
797 ++n;
798 }
799 }
800 if (*flag){
801 for (k=1;k<ktmax;++k){
802 for(j=0;j<jmax+2;++j){
803 for(i=0;i<imax+2;++i){
804 if (gridmap[j*(imax+2) + i]!=-1){
805 array_for_output[n*number_of_properties + 0]=height[k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
806 array_for_output[n*number_of_properties+1]=velocity[0][k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
807 array_for_output[n*number_of_properties+2]=velocity[1][k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
808 array_for_output[n*number_of_properties+3]=velocity[2][k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
809 array_for_output[n*number_of_properties+4]=drainage[k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
810 array_for_output[n*number_of_properties+5]=melt[k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
811 array_for_output[n*number_of_properties+6]=ice_land_sea[k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
812 array_for_output[n*number_of_properties+7]=type_of_base[k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
813 array_for_output[n*number_of_properties+8]=temperature[k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
814 array_for_output[n*number_of_properties+9]=age[k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
815 array_for_output[n*number_of_properties+10]=flux[0][k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
816 array_for_output[n*number_of_properties+11]=flux[1][k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
817 ++n;
818 }
819 }
820 }
821 }
822 kcmin=0;
823 }else{
824 kcmin=1;
825 }
826 for (k=kcmin;k<kcmax+1;++k){
827 for(j=0;j<jmax+2;++j){
828 for(i=0;i<imax+2;++i){
829 if (gridmap[j*(imax+2) + i]!=-1){
830 array_for_output[n*number_of_properties + 0]=height[nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
831 array_for_output[n*number_of_properties+1]=velocity[0][nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
832 array_for_output[n*number_of_properties+2]=velocity[1][nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
833 array_for_output[n*number_of_properties+3]=velocity[2][nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
834 array_for_output[n*number_of_properties+4]=drainage[nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
835 array_for_output[n*number_of_properties+5]=melt[nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
836 array_for_output[n*number_of_properties+6]=ice_land_sea[nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
837 array_for_output[n*number_of_properties+7]=type_of_base[nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
838 array_for_output[n*number_of_properties+8]=temperature[nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
839 array_for_output[n*number_of_properties+9]=age[nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
840 array_for_output[n*number_of_properties+10]=flux[0][nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
841 array_for_output[n*number_of_properties+11]=flux[1][nodes_in_temp+k*nodes_in_layer_of_staggered_grid+j*(imax+2)+i];
842 ++n;
843 }
844 }
845 }
846 }
847 if (n!=number_of_written_nodes){
848 printf("ERROR: Number of nodes for data %d does not match calculated %d\n", n, number_of_written_nodes);
849 }else{
850 /* write binary data on file */
851 sprintf(data_filename,"%5s%2s%s",runname, ergnum, suffix);
852 printf("| Writing on file %s\n", data_filename);
853 if((ptFile=fopen(data_filename, "w"))==NULL){
854 printf("ERROR: Could not open Elmer timeslice %s file for writing!\n", data_filename);
855 }else{
856 fprintf(ptFile, "0 0 %d %d\n", number_of_written_nodes, number_of_properties);
857 auxiliary = fwrite((void *) array_for_output, (size_t) number_of_properties*sizeof(float), (size_t) number_of_written_nodes, ptFile);
858 if (auxiliary!=number_of_written_nodes){
859 printf("ERROR: Only %d of %d items written on file\n",auxiliary,number_of_nodes);
860 printf("ERROR: Not able to write binary data on time-slice file %s\n", data_filename);
861 }else{
862 printf("| Succeeded in writing file\n");
863 printf("---------------------------------------------------------------\n");
864 }
865 }
866 }
867 fclose(ptFile);
868 free(array_for_interpolation);
869 free(float_maske);
870 free(float_cts);
871 free(iced);
872 free(gridmap);
873 free(array_for_output);
874 return;
875 }
876 /******************************************************
877 Output of SICOPOLIS grid in ELMER Solver format
878 *******************************************************/
FC_FUNC(pregrid,PREGRID)879 void STDCALLBULL FC_FUNC(pregrid,PREGRID) (float *xi, /* unscaled x coordinate index i: 0,imax */
880 float *eta, /* unscaled y coordinate index j from 0,jmax */
881 float *z_c, /* unscaled z coordinate index i: 0,imax, j: 0,jmax, kc: 0,kcmax */
882 float *z_t, /* unscaled y coordinate index i: 0,imax, j: 0,jmax, kt: 0,kcmax */
883 int *imax_in, /* grid steps in xi-direction */
884 int *jmax_in, /* grid steps in eta-direction */
885 int *kcmax_in, /* grid steps in z-direction in cold ice layer */
886 int *ktmax_in, /* grid steps in z-direction in temperate ice layer */
887 FC_CHAR_PTR(runname,runname_l), /*name of run*/
888 FC_CHAR_PTR(ergnum,ergnum_l), /*number of file*/
889 int *maske, /*mask of vertex type */
890 float *deltaX, /* stepsize of grid */
891 int *flag)
892 {
893 register int i, j, k, n, element_nr, boundary_nr,kn;
894 int number_of_layers[2], number_of_bulk_elements, number_of_boundary_elements, elements_in_one_layer, number_of_iced_nodes, number_of_nodes[2], nodes_of_element[8], *nodes_of_side_element, *parent_of_side_element, number_of_iced_collums, nodes_in_one_layer, number_of_ice_boundaries, *iced, *boundary, *gridmap, auxiliary, idnr, min_max_i_j[2][2];
895 int imax, jmax, kcmax, ktmax;
896 float *staggered_grid[2], min_max_xyz[2][3];
897 float actual_scaled_coord[3];
898 char groupid[4], filename[80], yes_no, *suffix=".ep";
899 FILE *ptFile;
900
901 /* constants */
902 imax= *imax_in;
903 jmax= *jmax_in;
904 kcmax= *kcmax_in;
905 ktmax= *ktmax_in;
906 nodes_in_one_layer = (imax+2)*(jmax+2);
907 elements_in_one_layer = (imax+1)*(jmax+1);
908 number_of_layers[0] = kcmax+1;
909 number_of_layers[1] = (*flag)*(ktmax);
910 number_of_nodes[0] = nodes_in_one_layer*number_of_layers[0];
911 number_of_nodes[1] = nodes_in_one_layer*number_of_layers[1];
912 /* print out little summary */
913 printf("---------------------------------------------------------------\n");
914 printf("| Output of SICOPOLIS Grid for ELMER Solver\n");
915 printf("---------------------------------------------------------------\n");
916 printf("| imax/jmax/kcmax/ktmax=%d/%d/%d/%d\n",imax, jmax, kcmax, ktmax);
917 printf("---------------------------------------------------------------\n");
918 printf("| nodes in original grid:\n");
919 printf("| cold layer: %d=%d * %d\n", (imax+1)*(jmax+1)*(kcmax+1), (imax+1)*(jmax+1), (kcmax+1));
920 printf("| temperate layer: %d=%d * %d\n", (imax+1)*(jmax+1)*(ktmax+1), (imax+1)*(jmax+1), (ktmax+1));
921 printf("| -------------\n");
922 printf("| %d=%d * %d\n", (imax+1)*(jmax+1)*(ktmax+2+kcmax), (imax+1)*(jmax+1), ktmax+2+kcmax);
923 printf("| x = %.1f -> %.1f , y = %.1f -> %.1f, dx=%.1f\n", xi[0], xi[imax], eta[0], eta[jmax], *deltaX);
924 printf("---------------------------------------------------------------\n");
925 if (*flag)
926 printf("| Output of temperate layer enabled\n");
927 else printf("| Output of temperate layer disabled\n");
928 printf("---------------------------------------------------------------\n");
929 printf("| nodes in full staggered grid:\n");
930 printf("| cold layer: %d=%d*%d\n", number_of_nodes[0], nodes_in_one_layer, number_of_layers[0]);
931 printf("| temperate layer: %d=%d*%d\n", number_of_nodes[1], nodes_in_one_layer, number_of_layers[1]);
932 printf("| -------------\n");
933 printf("| %d\n", number_of_nodes[1] + number_of_nodes[0]);
934 printf("---------------------------------------------------------------\n");
935 /* allocate memory */
936 staggered_grid[0] = (float *) malloc((size_t) (number_of_nodes[0]*3*sizeof(float)));
937 if (staggered_grid[0] == NULL){
938 printf("ERROR in allocating memory for staggered grid of cold ice layer\n");
939 return;
940 }
941 staggered_grid[1] = (float *) malloc((size_t) (number_of_nodes[1]*3*sizeof(float)));
942 if (staggered_grid[1] == NULL){
943 printf("ERROR in allocating memory for staggered grid of temperate ice layer\n");
944 free(staggered_grid[0]);
945 return;
946 }
947 iced = (int *) malloc((size_t) ((size_t) (imax+1)*(jmax+1)*sizeof(int)));
948 if (iced == NULL){
949 printf("ERROR in allocating memory for glaciation information\n");
950 free(staggered_grid[0]);free(staggered_grid[1]);
951 return;
952 }
953 /* get staggered grid(s) */
954 auxiliary = get_staggered_grid(xi,eta,z_c,imax,jmax,kcmax,deltaX,staggered_grid[0]);
955 if (auxiliary!=number_of_nodes[0]){
956 printf("WARNING: number of written %d gridpoints for cold layer does not match calculated %d\n",auxiliary, number_of_nodes[0]);
957 printf("ERROR: Interpolation of staggered grid for cold layer failed!\n");
958 free(staggered_grid[0]);free(staggered_grid[1]);
959 return;
960 }else
961 printf("| succeeded in interpolating staggered grid for cold layer\n");
962 if (*flag){
963 auxiliary = get_staggered_grid(xi,eta,z_t,imax,jmax,(ktmax-1),deltaX,staggered_grid[1]);
964 if(auxiliary!=number_of_nodes[1]){
965 printf("WARNING: number of written %d gridpoints for cold layer does not match calculated %d\n",auxiliary, number_of_nodes[1]);
966 printf("ERROR: Interpolation of staggered grid for temperate layer failed!\n");
967 free(staggered_grid[0]);free(staggered_grid[1]);
968 return;
969 }else
970 printf("| succeeded in interpolating staggered grid for temperate layer\n");
971 }else
972 printf("| no staggered grid for temperate layer interpolated\n");
973 printf("---------------------------------------------------------------\n");
974 printf("| staggered grid:\n");
975 printf("| x = %.1f -> %.1f , y = %.1f -> %.1f\n",staggered_grid[0][0], staggered_grid[0][(imax+1)*3], staggered_grid[0][1], staggered_grid[0][(jmax+1)*(imax+2)*3+1]);
976 printf("---------------------------------------------------------------\n");
977 /* get glaciation info */
978 number_of_iced_collums=get_glaciation_info(imax,jmax,iced,maske);
979 if (number_of_iced_collums<1){
980 printf("| no glaciation\n");
981 /* printf("ERROR in calculation of glaciation\n"); */
982 /* free(staggered_grid[0]);free(staggered_grid[1]);free(iced); */
983 /* return; */
984 }else{
985 printf("| number of iced colums: %d out of %d (%3.2f %%)\n", number_of_iced_collums, (imax+1)*(jmax+1), ((float) number_of_iced_collums)*100.0/((float) (imax+1)*(jmax+1)));
986 }
987 gridmap = (int *) malloc((size_t) nodes_in_one_layer*sizeof(int));
988 if (gridmap==NULL){
989 printf("ERROR in allocation of memory for staggered grid glaciation information\n");
990 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);
991 return;
992 }
993 number_of_iced_nodes = get_grid_map(imax,jmax,iced,gridmap);
994 /* if (number_of_iced_nodes<1){ */
995 /* printf("ERROR in calculation of glaciation information for staggered grid\n"); */
996 /* free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(gridmap); */
997 /* return; */
998 /* } */
999 printf("| number of iced nodes in one layer: %d out of %d\n",number_of_iced_nodes, nodes_in_one_layer);
1000 boundary = (int *) malloc((size_t) number_of_iced_collums*4*sizeof(int));
1001 if (boundary==NULL){
1002 printf("ERROR in allocation of memory for ice-boundary information\n");
1003 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(gridmap);
1004 return;
1005 }
1006 number_of_ice_boundaries = get_glaciation_boundary_info(imax,jmax,iced,boundary);
1007 /* if (number_of_ice_boundaries<1){ */
1008 /* printf("ERROR in calculation of boundaries\n"); */
1009 /* free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap); */
1010 /* return; */
1011 /* } */
1012 printf("| number of boundary elements lines: %d (%f km)\n", number_of_ice_boundaries, ((float) number_of_ice_boundaries)*(*deltaX));
1013 /* calculate constants for output mesh*/
1014 if (*flag){
1015 number_of_nodes[0] = number_of_iced_nodes*(kcmax+1);
1016 number_of_nodes[1] = number_of_iced_nodes*(ktmax);
1017 }else{
1018 number_of_nodes[0] = number_of_iced_nodes*(kcmax+1);
1019 number_of_nodes[1] = 0;
1020 }
1021 number_of_bulk_elements = (kcmax+(*flag)*ktmax)*number_of_iced_collums;
1022 number_of_boundary_elements = number_of_ice_boundaries*(kcmax+(*flag)*ktmax) + (2+(*flag))*number_of_iced_collums;
1023 nodes_of_side_element = (int *) malloc((size_t) number_of_ice_boundaries*4*(kcmax+(*flag)*ktmax)*sizeof(int));
1024 if (nodes_of_side_element == NULL){
1025 printf("ERROR in allocation of side element array\n");
1026 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);
1027 return;
1028 }
1029 parent_of_side_element = (int *) malloc((size_t) number_of_ice_boundaries*(kcmax+(*flag)*ktmax)*sizeof(int));
1030 if (parent_of_side_element == NULL){
1031 printf("ERROR in allocation of side element array\n");
1032 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);
1033 return;
1034 }
1035 printf("---------------------------------------------------------------\n");
1036 printf("| number of nodes:\n");
1037 printf("| in cold layer %d\n",number_of_nodes[0]);
1038 printf("| in temperate layer %d\n",number_of_nodes[1]);
1039 printf("| -----------\n");
1040 printf("| %d\n",number_of_nodes[0]+number_of_nodes[1]);
1041 printf("---------------------------------------------------------------\n");
1042 printf("| number of elements:\n");
1043 printf("| in cold layer volume %d\n",kcmax*number_of_iced_collums);
1044 printf("| in temperate layer volume %d\n",(*flag)*ktmax*number_of_iced_collums);
1045 printf("| -----------\n");
1046 printf("| %d\n",number_of_bulk_elements);
1047 printf("| -----------\n");
1048 printf("| \n");
1049 printf("| at base %d\n",number_of_iced_collums);
1050 printf("| at cts %d\n",number_of_iced_collums);
1051 printf("| at free surface %d\n",number_of_iced_collums);
1052 printf("| on boundary of cold layer %d\n",number_of_ice_boundaries*(kcmax));
1053 printf("| on boundary of temperate layer %d\n",(*flag)*number_of_ice_boundaries*(ktmax));
1054 printf("| -----------\n");
1055 printf("| %d\n",number_of_boundary_elements);
1056 printf("---------------------------------------------------------------\n");
1057 /* writing header file */
1058 sprintf(filename,"mesh.header");
1059 printf("| Writing mesh header file %s.\n",filename);
1060 if((ptFile=fopen(filename, "w"))==NULL){
1061 printf("\a Could not open Elmer mesh-file file for writing!\n");
1062 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);free(parent_of_side_element);
1063 return;
1064 }
1065 fprintf(ptFile,"%d %d %d\n",number_of_nodes[0]+number_of_nodes[1], number_of_bulk_elements,number_of_boundary_elements);
1066 fprintf(ptFile,"2\n");
1067 fprintf(ptFile,"808 %d\n", number_of_bulk_elements);
1068 fprintf(ptFile,"404 %d\n", number_of_boundary_elements);
1069 printf("| succeeded in writting mesh header file %s.\n",filename);
1070 printf("---------------------------------------------------------------\n");
1071 fclose(ptFile);
1072 /* check min/max values of grid */
1073 min_max_i_j[0][0]=imax;
1074 min_max_i_j[1][0]=0;
1075 min_max_i_j[0][1]=jmax;
1076 min_max_i_j[1][1]=0;
1077 for(j=0;j<jmax+1;++j){
1078 for (i=0;i<imax+1;++i){
1079 if (iced[j*(imax+1)+i]!=-1){
1080 if (i < min_max_i_j[0][0]) min_max_i_j[0][0]=i;
1081 if (i > min_max_i_j[1][0]) min_max_i_j[1][0]=i;
1082 if (j < min_max_i_j[0][1]) min_max_i_j[0][1]=j;
1083 if (j > min_max_i_j[1][1]) min_max_i_j[1][1]=j;
1084 }
1085 }
1086 }
1087 printf("| glaciation information:\n");
1088 printf("| original grid:\n");
1089 printf("| i = %d -> %d of %d \n", min_max_i_j[0][0], min_max_i_j[1][0], imax);
1090 printf("| j = %d -> %d of %d \n", min_max_i_j[0][1], min_max_i_j[1][1], jmax);
1091 min_max_i_j[0][0]=imax+1;
1092 min_max_i_j[1][0]=0;
1093 min_max_i_j[0][1]=jmax+1;
1094 min_max_i_j[1][1]=0;
1095 for (i=0;i<(imax+1)*(jmax+1);++i){
1096 if (gridmap[i]!=-1){
1097 for (n=0;n<3;++n)
1098 min_max_xyz[0][n]=staggered_grid[0][i*3 + n];
1099 break;
1100 }
1101 }
1102 for(j=0;j<jmax+2;++j){
1103 for (i=0;i<imax+2;++i){
1104 if (gridmap[j*(imax+2) + i]!=-1){
1105 if (i < min_max_i_j[0][0]) min_max_i_j[0][0]=i;
1106 if (i > min_max_i_j[1][0]) min_max_i_j[1][0]=i;
1107 if (j < min_max_i_j[0][1]) min_max_i_j[0][1]=j;
1108 if (j > min_max_i_j[1][1]) min_max_i_j[1][1]=j;
1109 for (n=0;n<2;++n){
1110 if (min_max_xyz[0][n] > staggered_grid[0][(j*(imax+2) + i)*3 + n]) min_max_xyz[0][n] = staggered_grid[0][(j*(imax+2) + i)*3 + n];
1111 if (min_max_xyz[1][n] < staggered_grid[0][(j*(imax+2) + i)*3 + n]) min_max_xyz[1][n] = staggered_grid[0][(j*(imax+2) + i)*3 + n];
1112 }
1113 if (*flag){
1114 if (min_max_xyz[0][2] > staggered_grid[0][(j*(imax+2) + i)*3 + 2]) min_max_xyz[0][2] = staggered_grid[1][(j*(imax+2) + i)*3 + 2];
1115 }
1116 for (k=0;k<(kcmax+1);++k){
1117 if (min_max_xyz[0][2] > staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 2]) min_max_xyz[0][2] = staggered_grid[0][(j*(imax+2) + i)*3 + 2];
1118 if (min_max_xyz[1][2] < staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 2]) min_max_xyz[1][2] = staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 2];
1119 }
1120 }
1121 }
1122 }
1123 printf("| staggered grid:\n");
1124 printf("| i = %d -> %d of %d \n", min_max_i_j[0][0], min_max_i_j[1][0], imax);
1125 printf("| j = %d -> %d of %d \n", min_max_i_j[0][1], min_max_i_j[1][1], jmax);
1126 printf("| x = %.8f -> %.8f\n",min_max_xyz[0][0], min_max_xyz[1][0]);
1127 printf("| y = %.8f -> %.8f\n",min_max_xyz[0][1], min_max_xyz[1][1]);
1128 printf("| z = %.8f -> %.8f\n",min_max_xyz[0][2], min_max_xyz[1][2]);
1129 printf("---------------------------------------------------------------\n");
1130
1131 /* writing nodes file */
1132 sprintf(filename,"mesh.nodes");
1133 printf("| Writing node-file %s.\n",filename);
1134 if((ptFile=fopen(filename, "w"))==NULL){
1135 printf("\a Could not open Elmer mesh-file file for writing!\n");
1136 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);free(parent_of_side_element);
1137 return;
1138 }
1139 n=0;
1140 if (*flag){
1141 for (k=0;k<ktmax;++k){
1142 for(j=0;j<jmax+2;++j){
1143 for (i=0;i<imax+2;++i){
1144 if (gridmap[j*(imax+2) + i]!=-1){
1145 fprintf(ptFile, "%d -1 %.8f %.8f %.8f\n", n+1, staggered_grid[1][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 0], staggered_grid[1][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 1], staggered_grid[1][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 2]);
1146 ++n;
1147 }
1148 }
1149 }
1150 }
1151 }
1152 for (k=0;k<kcmax+1;++k){
1153 for(j=0;j<jmax+2;++j){
1154 for (i=0;i<imax+2;++i){
1155 if (gridmap[j*(imax+2) + i]!=-1){
1156 fprintf(ptFile, "%d -1 %.8f %.8f %.8f\n", n+1, staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 0], staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 1], staggered_grid[0][(k*nodes_in_one_layer + j*(imax+2) + i)*3 + 2]);
1157 ++n;
1158 }
1159 }
1160 }
1161 }
1162 if (n!=number_of_nodes[0]+number_of_nodes[1]){
1163 printf("ERROR: number of written nodes %d does not match calculated %d\n",n, number_of_nodes[0]+number_of_nodes[1]);
1164 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);free(parent_of_side_element);
1165 fclose(ptFile);
1166 return;
1167 }
1168 printf("| succeeded in writting node file %s.\n",filename);
1169 printf("---------------------------------------------------------------\n");
1170 /* writing element file */
1171 sprintf(filename,"mesh.elements");
1172 printf("| Writing bulk element file %s.\n",filename);
1173 if((ptFile=fopen(filename, "w"))==NULL){
1174 printf("\a Could not open Elmer mesh-file file for writing!\n");
1175 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);free(parent_of_side_element);
1176 return;
1177 }
1178 n=0;element_nr=0;boundary_nr=0; idnr=1;
1179 if (*flag){
1180 for (k=0;k<ktmax;++k){
1181 for (j=0;j<jmax+1;++j){
1182 for (i=0;i<imax+1;++i){
1183 if (iced[j*(imax+1)+i]!=-1){
1184 for (n=0; n<2; n++){/* lower level: n=0; upper level n=1; each counterclkws */
1185 nodes_of_element[n*4] = (k+n)*number_of_iced_nodes + gridmap[j*(imax+2)+i] + 1;
1186 nodes_of_element[n*4+1] = (k+n)*number_of_iced_nodes + gridmap[j*(imax+2)+i+1] + 1;
1187 nodes_of_element[n*4+2] = (k+n)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i+1] + 1;
1188 nodes_of_element[n*4+3] = (k+n)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i] + 1;
1189 }
1190 fprintf(ptFile,"%d %d 808 %d %d %d %d %d %d %d %d\n", element_nr+1, idnr,
1191 nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3],
1192 nodes_of_element[4], nodes_of_element[5], nodes_of_element[6], nodes_of_element[7]);
1193 if (element_nr != k+number_of_iced_collums+ iced[j*(imax+1)+i]) printf("element_nr=%d does not match %d\n", element_nr, k+number_of_iced_collums+ iced[j*(imax+1)+i]);
1194 ++element_nr;
1195 if (boundary[iced[(j*(imax+1)+i)]*4+0]){/* east */
1196 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[2];
1197 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[1];
1198 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[5];
1199 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[6];
1200 parent_of_side_element[boundary_nr]=element_nr;
1201 ++boundary_nr;
1202 }
1203 if (boundary[iced[(j*(imax+1)+i)]*4+1]){/* north */
1204 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
1205 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[2];
1206 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[6];
1207 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
1208 parent_of_side_element[boundary_nr]=element_nr;
1209 ++boundary_nr;
1210 }
1211 if (boundary[iced[(j*(imax+1)+i)]*4+2]){/* west */
1212 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
1213 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
1214 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
1215 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
1216 parent_of_side_element[boundary_nr]=element_nr;
1217 ++boundary_nr;
1218 }
1219 if (boundary[iced[(j*(imax+1)+i)]*4+3]){/* south */
1220 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[1];
1221 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
1222 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
1223 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[5];
1224 parent_of_side_element[boundary_nr]=element_nr;
1225 ++boundary_nr;
1226 }
1227 }
1228 }
1229 }
1230 }
1231 idnr=2;
1232 }
1233 for (k=0;k<kcmax;++k){
1234 for (j=0;j<jmax+1;++j){
1235 for (i=0;i<imax+1;++i){
1236 if (iced[j*(imax+1)+i]!=-1){
1237 for (n=0; n<2; n++){/* lower level: n=0; upper level n=1; each counterclkws */
1238 nodes_of_element[n*4] = number_of_nodes[1] + (k+n)*number_of_iced_nodes + gridmap[j*(imax+2)+i] + 1;
1239 nodes_of_element[n*4+1] = number_of_nodes[1] + (k+n)*number_of_iced_nodes + gridmap[j*(imax+2)+i+1] + 1;
1240 nodes_of_element[n*4+2] = number_of_nodes[1] + (k+n)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i+1] + 1;
1241 nodes_of_element[n*4+3] = number_of_nodes[1] + (k+n)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i] + 1;
1242 }
1243 fprintf(ptFile,"%d %d 808 %d %d %d %d %d %d %d %d\n", element_nr+1, idnr,
1244 nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3],
1245 nodes_of_element[4], nodes_of_element[5], nodes_of_element[6], nodes_of_element[7]);
1246 if (element_nr != k*number_of_iced_collums+ iced[j*(imax+1)+i]) printf("element_nr=%d does not match %d\n", element_nr, k*number_of_iced_collums+ iced[j*(imax+1)+i]);
1247 ++element_nr;
1248 if (boundary[iced[(j*(imax+1)+i)]*4+0]){/* east */
1249 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[2];
1250 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[1];
1251 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[5];
1252 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[6];
1253 parent_of_side_element[boundary_nr]=element_nr;
1254 ++boundary_nr;
1255 }
1256 if (boundary[iced[(j*(imax+1)+i)]*4+1]){/* north */
1257 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
1258 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[2];
1259 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[6];
1260 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
1261 parent_of_side_element[boundary_nr]=element_nr;
1262 ++boundary_nr;
1263 }
1264 if (boundary[iced[(j*(imax+1)+i)]*4+2]){/* west */
1265 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[3];
1266 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
1267 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
1268 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[7];
1269 parent_of_side_element[boundary_nr]=element_nr;
1270 ++boundary_nr;
1271 }
1272 if (boundary[iced[(j*(imax+1)+i)]*4+3]){/* south */
1273 nodes_of_side_element[boundary_nr*4+0]=nodes_of_element[1];
1274 nodes_of_side_element[boundary_nr*4+1]=nodes_of_element[0];
1275 nodes_of_side_element[boundary_nr*4+2]=nodes_of_element[4];
1276 nodes_of_side_element[boundary_nr*4+3]=nodes_of_element[5];
1277 parent_of_side_element[boundary_nr]=element_nr;
1278 ++boundary_nr;
1279 }
1280 }
1281 }
1282 }
1283 }
1284 if (number_of_bulk_elements!=element_nr){
1285 printf("ERROR: Number of written bulk elements %d does not match calculated %d\n", number_of_bulk_elements, element_nr);
1286 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);free(parent_of_side_element);
1287 return;
1288 }
1289 printf("| succeeded in writting bulk element file %s.\n",filename);
1290 printf("---------------------------------------------------------------\n");
1291 fclose(ptFile);
1292 /* writing boundary element file */
1293 sprintf(filename,"mesh.boundary");
1294 printf("| Writing boundary element file %s.\n",filename);
1295 if((ptFile=fopen(filename, "w"))==NULL){
1296 printf("\a Could not open Elmer boundary element file for writing!\n");
1297 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);free(parent_of_side_element);
1298 return;
1299 }
1300 boundary_nr=0;
1301 /* base */
1302 idnr=1;
1303 for (j=0;j<jmax+1;++j){
1304 for (i=0;i<imax+1;++i){
1305 if (iced[j*(imax+1)+i]!=-1){
1306 nodes_of_element[0]= gridmap[j*(imax+2)+i]+1;
1307 nodes_of_element[1]= gridmap[j*(imax+2)+i+1]+1;
1308 nodes_of_element[2]= gridmap[(j+1)*(imax+2)+i+1]+1;
1309 nodes_of_element[3]= gridmap[(j+1)*(imax+2)+i]+1;
1310 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n", boundary_nr+1, idnr, iced[j*(imax+1)+i]+1, nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3]);
1311 ++boundary_nr;
1312 }
1313 }
1314 }
1315 /* cts */
1316 if (*flag){
1317 ++idnr;
1318 for (j=0;j<jmax+1;++j){
1319 for (i=0;i<imax+1;++i){
1320 if (iced[j*(imax+1)+i]!=-1){
1321 nodes_of_element[0] = number_of_nodes[1] + gridmap[j*(imax+2)+i]+1;
1322 nodes_of_element[1] = number_of_nodes[1] + gridmap[j*(imax+2)+i+1]+1;
1323 nodes_of_element[2] = number_of_nodes[1] + gridmap[(j+1)*(imax+2)+i+1]+1;
1324 nodes_of_element[3] = number_of_nodes[1] + gridmap[(j+1)*(imax+2)+i]+1;
1325 fprintf(ptFile,"%d %d %d %d 404 %d %d %d %d\n", boundary_nr+1, idnr, (ktmax-2)*number_of_iced_collums+iced[j*(imax+1)+i]+1, (ktmax)*number_of_iced_collums+iced[j*(imax+1)+i]+1, nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3]);
1326 ++boundary_nr;
1327 }
1328 }
1329 }
1330 }
1331 /* free surface */
1332 ++idnr;
1333 for (j=0;j<jmax+1;++j){
1334 for (i=0;i<imax+1;++i){
1335 if (iced[j*(imax+1)+i]!=-1){
1336 nodes_of_element[0] = number_of_nodes[1] + (kcmax)*number_of_iced_nodes + gridmap[j*(imax+2)+i]+1;
1337 nodes_of_element[1] = number_of_nodes[1] + (kcmax)*number_of_iced_nodes + gridmap[j*(imax+2)+i+1]+1;
1338 nodes_of_element[2] = number_of_nodes[1] + (kcmax)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i+1]+1;
1339 nodes_of_element[3] = number_of_nodes[1] + (kcmax)*number_of_iced_nodes + gridmap[(j+1)*(imax+2)+i]+1;
1340 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n", boundary_nr+1, idnr, ((ktmax)*(*flag)+kcmax-1)*number_of_iced_collums+iced[j*(imax+1)+i]+1, nodes_of_element[0], nodes_of_element[1], nodes_of_element[2], nodes_of_element[3]);
1341 ++boundary_nr;
1342 }
1343 }
1344 }
1345 /* sides */
1346 if (*flag){
1347 ++idnr;
1348 for (n=0;n<number_of_ice_boundaries*ktmax;++n){
1349 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n", boundary_nr+1, parent_of_side_element[n], idnr, nodes_of_side_element[n*4+0], nodes_of_side_element[n*4+1], nodes_of_side_element[n*4+2], nodes_of_side_element[n*4+3]);
1350 ++boundary_nr;
1351 }
1352 ++idnr;
1353 for (;n<number_of_ice_boundaries*(kcmax+ktmax);++n){
1354 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n", boundary_nr+1, idnr, parent_of_side_element[n], nodes_of_side_element[n*4+0], nodes_of_side_element[n*4+1], nodes_of_side_element[n*4+2], nodes_of_side_element[n*4+3]);
1355 ++boundary_nr;
1356 }
1357 }else{
1358 ++idnr;
1359 for (n=0;n<number_of_ice_boundaries*kcmax;++n){
1360 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n", boundary_nr+1, idnr, parent_of_side_element[n], nodes_of_side_element[n*4+0], nodes_of_side_element[n*4+1], nodes_of_side_element[n*4+2], nodes_of_side_element[n*4+3]);
1361 ++boundary_nr;
1362 }
1363 }
1364 if (boundary_nr!=number_of_boundary_elements){
1365 printf("ERROR: Number of written boundary elements %d does not match calculated %d.\n", boundary_nr, number_of_boundary_elements);
1366 }else{
1367 printf("| succeeded in writting boundary element file %s.\n",filename);
1368 printf("---------------------------------------------------------------\n");
1369 }
1370 fclose(ptFile);
1371 free(staggered_grid[0]);free(staggered_grid[1]);free(iced);free(boundary);free(gridmap);free(nodes_of_side_element);free(parent_of_side_element);
1372 return;
1373 }
1374
FC_FUNC(asciidata,ASCIIDATA)1375 void STDCALLBULL FC_FUNC(asciidata,ASCIIDATA) (float *xi, /* unscaled x coordinate index i: 0,imax */
1376 float *eta, /* unscaled y coordinate index j from 0,jmax */
1377 int *imax_in, /* grid steps in xi-direction */
1378 int *jmax_in, /* grid steps in eta-direction */
1379 int *kcmax_in, /* grid steps in z-direction in cold ice layer */
1380 int *ktmax_in, /* grid steps in z-direction in temperate ice layer */
1381 float *z_c, /* z-coordinate in cold layer for given i,j-position in plane and kc in vertical*/
1382 float *z_t, /* z-coordinate in temperate layer for given i,j-position in plane and kc in vertical*/
1383 float *vx_c, /* velocity component in x for cold region for given i,j-position in plane and kc in vertical */
1384 float *vy_c, /* velocity component in y for cold region for given i,j-position in plane and kc in vertical */
1385 float *vz_c, /* velocity component in z for cold region for given i,j-position in plane and kc in vertical */
1386 float *age_c, /* age in cold region for given i,j-position in plane and kc in vertical */
1387 float *temp_c, /* temperature in cold region for given i,j-position in plane and kt in vertical */
1388 float *vx_t, /* velocity component in x for temperated region for given i,j-position in plane and kt in vertical */
1389 float *vy_t, /* velocity component in y for temperated region for given i,j-position in plane and kt in vertical */
1390 float *vz_t, /* velocity component in z for temperated region for given i,j-position in plane and kt in vertical */
1391 float *temp_t_m, /* melting temperature in temperate ice region for given i,j-position in plane and kt in vertical */
1392 float *age_t, /* age in temperate region for given i,j-position in plane and kc in vertical */
1393 float *omega_t, /* H2O mass fraction in temperated region for given i,j-position in plane and kt in vertical */
1394 float *Q_bm, /* production rate of melting-water at bottom for given i,j-position in plane */
1395 float *Q_tld, /* water drainage rate from the temperated region for given i,j-position in plane */
1396 float *am_perp, /* ice volume flux through CTS for given i,j-position in plane */
1397 float *qx, /* mass-flux in x-direction for given i,j-position in plane */
1398 float *qy, /* mass-flux in y-direction for given i,j-position in plane */
1399 int *n_cts, /* polythermal conditions for given i,j-position at base (-1 cold ice; 0 temp. ice base with cold ice above; 1 temp. ice base with temperate ice above; */
1400 int *maske, /* glaciation information for given i,j-position at base (glaciated=0, ice-free=1, 2=sea-point, 3=floating ice) */
1401 FC_CHAR_PTR(runname,p1), /*name of run*/
1402 FC_CHAR_PTR(ergnum,p2), /*number of file*/
1403 int *flag)
1404 {
1405 register int i, j, k, n, current_index;
1406 int imax, jmax, kcmax, ktmax, kcmin, offset, nodes_in_temp, nodes_in_cold, elements_in_layer, number_of_iced_nodes, number_of_written_nodes, *gridmap, ok;
1407 int nodes_in_one_layer, number_of_nodes, nodes_in_layer_of_staggered_grid, number_of_iced_collums, *iced, auxiliary, number_of_properties, outputflags[12], failed=0;
1408 char filename[80], inputstring[40], dummy_string[39], user_name[10],*suffix=".asc";
1409 struct stat buf;
1410 time_t how_late_is_it;
1411 FILE *ptFile;
1412 /* constants */
1413 imax=*imax_in;
1414 jmax=*jmax_in;
1415 kcmax=*kcmax_in;
1416 ktmax=*ktmax_in;
1417 elements_in_layer=(imax+1)*(jmax+1);
1418 printf("---------------------------------------------------------------\n");
1419 printf("| Output of SICOPOLIS Data in ASCII-format - Hi Olli :-)\n");
1420 printf("---------------------------------------------------------------\n");
1421 printf("| imax/jmax/kcmax/ktmax=%d/%d/%d/%d\n",imax, jmax, kcmax, ktmax);
1422 printf("---------------------------------------------------------------\n");
1423 printf("| nodes in grid:\n");
1424 printf("| cold layer: %d=%d * %d\n", (imax+1)*(jmax+1)*(kcmax+1), (imax+1)*(jmax+1), (kcmax+1));
1425 printf("| temperate layer: %d=%d * %d\n", (imax+1)*(jmax+1)*(ktmax+1), (imax+1)*(jmax+1), (ktmax+1));
1426 printf("| -------------\n");
1427 printf("| %d=%d * %d\n", (imax+1)*(jmax+1)*(ktmax+2+kcmax), (imax+1)*(jmax+1), ktmax+2+kcmax);
1428 printf("---------------------------------------------------------------\n");
1429 printf("| Trying to load preferences from .sico2elmer.rc\n");
1430 /* read preference file */
1431 if (stat(".sico2elmer.rc" , &buf)!=0){
1432 if (errno == ENOENT){
1433 printf("WARNING: The preference-file .sico2elmer.rc does not exist. Writing all output-components by default\n");
1434 failed=1;
1435 }
1436 else if (errno == EACCES)
1437 printf("WARNING: No permission to read file .sico2elmer.rc. Writing all output-components by default\n");
1438 else printf("WARNING: Error occured during reading preference-file .sico2elmer.rc. Writing all output-components by default\n");
1439 failed=1;
1440 }else{
1441 if((ptFile=fopen(".sico2elmer.rc", "r"))==NULL){
1442 printf("WARNING: Error while opening file .sico2elmer.rc . Writing all output-components\n");
1443 failed=1;
1444 }else{
1445 /* read in comment lines */
1446 for (n=0;n<7;++n){
1447 if (fgets(inputstring, 40, ptFile)==NULL){
1448 printf("WARNING: Error while reading .sico2elmer.rc . Writing all output-components\n");
1449 failed=1;
1450 }
1451 }
1452 /* read: flag for output of temperate layer */
1453 if (fgets(inputstring, 40, ptFile)==NULL){
1454 printf("WARNING: Error while reading .sico2elmer.rc . Writing all output-components\n");
1455 failed=1;
1456 }else{
1457 if (sscanf(inputstring, "%1d%s",&outputflags[0], dummy_string)<0){
1458 printf("WARNING: Error while copying inputstring no 0 from .sico2elmer.rc\n");
1459 outputflags[0]=1;
1460 }
1461 }
1462 /* read: flag for output on not iced vertices*/
1463 if (fgets(inputstring, 40, ptFile)==NULL){
1464 printf("WARNING: Error while reading .sico2elmer.rc . Writing all output-components\n");
1465 failed=1;
1466 }else{
1467 if (sscanf(inputstring, "%1d%s",&outputflags[1], dummy_string)<0){
1468 printf("WARNING: Error while copying inputstring no 0 from .sico2elmer.rc\n");
1469 outputflags[1]=0;
1470 }
1471 }
1472 if (fgets(inputstring, 40, ptFile)==NULL){
1473 printf("WARNING: Error while reading .sico2elmer.rc. Writing all output-components\n");
1474 failed=1;
1475 }
1476 /* read: 3d variable flags */
1477 for (n=0;n<4;++n){
1478 if (fgets(inputstring, 40, ptFile)==NULL){
1479 printf("WARNING: Error while reading .sico2elmer.rc. Writing all output-components\n");
1480 failed=1;
1481 }else{
1482 if (sscanf(inputstring, "%1d%s",&outputflags[2+n], dummy_string)<0){
1483 printf("WARNING: Error while copying inputstring no %d from .sico2elmer.rc\n", n+1);
1484 outputflags[2+n]=1;
1485 }
1486 }
1487 }
1488 if (fgets(inputstring, 40, ptFile)==NULL){
1489 printf("WARNING: Error while reading .sico2elmer.rc. Writing all output-components\n");
1490 failed=1;
1491 }
1492 /* read: 2d variable flags */
1493 for (n=0;n<6;++n){
1494 if (fgets(inputstring, 40, ptFile)==NULL){
1495 printf("WARNING: Error while reading .sico2elmer.rc. Writing all output-components\n");
1496 failed=1;
1497 }else{
1498 if (sscanf(inputstring, "%1d%s",&outputflags[6+n], dummy_string)<0){
1499 printf("WARNING: Error while copying inputstring no %d from .sico2elmer.rc\n", n+1);
1500 outputflags[6+n]=1;
1501 }
1502 }
1503 }
1504 }
1505 fclose(ptFile);
1506 }
1507 if (failed){
1508 for (n=0;n<12;++n) outputflags[n]=1;
1509 printf("| Error occured during reading of file\n");
1510 printf("| Taking default values\n");
1511 }else printf("| Loading of preferences succeessful\n");
1512 printf("| Following values for output are set (1=yes,0=no):\n");
1513 printf("| \n");
1514 if (outputflags[0] && ((*flag)==0)){
1515 outputflags[0] = 0;
1516 printf("| Output temperate layer %d (reset, since no temperate layer has been loaded)\n", outputflags[0]);
1517 }else
1518 printf("| Output temperate layer %d\n", outputflags[0]);
1519 printf("| Output at not iced vertices %d\n", outputflags[1]);
1520 printf("| \n");
1521 printf("| Output 3d coordinates %d\n", outputflags[2]);
1522 printf("| Output velocities %d\n", outputflags[3]);
1523 printf("| Output temperature %d\n", outputflags[4]);
1524 printf("| Output age %d\n", outputflags[5]);
1525 printf("| \n");
1526 printf("| Output bedrock and ice depht %d\n", outputflags[6]);
1527 printf("| Output water drainage %d\n", outputflags[7]);
1528 printf("| Output mask for glaciation %d\n", outputflags[8]);
1529 printf("| Output type of base %d\n", outputflags[9]);
1530 printf("| Output flux %d\n", outputflags[10]);
1531 printf("| Output 2d coordinates %d\n", outputflags[11]);
1532 printf("---------------------------------------------------------------\n");
1533 /* open file for writing of 3d variables */
1534 sprintf(filename,"%5s%2s_3d%s", runname, ergnum, suffix);
1535 if((ptFile=fopen(filename, "w"))==NULL){
1536 printf("\a Could not open file for writing of 3d ASCII-data\n");
1537 return;
1538 }else{
1539 printf("| writing on 3d output file %s\n", filename);
1540 fprintf(ptFile,"# *********************************************************\n");
1541 fprintf(ptFile,"# ASCCI output file of run %s, file no. %s\n", runname, ergnum);
1542 if (time(&how_late_is_it)==-1){
1543 printf("WARNING: Could not evaluate current time\n");
1544 fprintf(ptFile,"# no specific date was able to be inquiered\n");
1545 }else{
1546 fprintf(ptFile,"# written on %s", ctime(&how_late_is_it));
1547 }
1548 /* replace with getpwuid(geteuid()) */
1549 #if defined(__APPLE__) || defined(MINGW32) || defined(WIN32) || defined(BSD)
1550 if ( 1 )
1551 #else
1552 if (cuserid(user_name)==NULL)
1553 #endif
1554 fprintf(ptFile,"# no user id was able to be inquiered\n");
1555 else
1556 fprintf(ptFile,"# by %10s\n", user_name);
1557 fprintf(ptFile,"# *********************************************************\n");
1558 fprintf(ptFile,"# 3d variables:\n");
1559 fprintf(ptFile,"# ");
1560 if (outputflags[1])
1561 fprintf(ptFile,"x y z ");
1562 if (outputflags[3])
1563 fprintf(ptFile,"v_x v_y v_z ");
1564 if (outputflags[4])
1565 fprintf(ptFile,"T ");
1566 if (outputflags[5])
1567 fprintf(ptFile,"Age ");
1568 if ((outputflags[2] + outputflags[3] + outputflags[4] + outputflags[5]) !=0){
1569 fprintf(ptFile,"\n# --------- temperate layer --------------------------------------------------------------------\n");
1570 if (outputflags[0]){
1571 for (k=0;k<ktmax;++k){
1572 for (j=0;j<jmax+1;++j){
1573 for (i=0;i<imax+1;++i){
1574 if( (outputflags[1]) && !((maske[j*(imax+1) + i]==0) || (maske[j*(imax+1) + i]==3)) ) continue;
1575 current_index = k*elements_in_layer + j*(imax+1) + i;
1576 if (outputflags[2])
1577 fprintf(ptFile,"% .4e % .4e % .4e ",xi[i], eta[j], z_t[current_index]);
1578 if (outputflags[3])
1579 fprintf(ptFile,"% .4e % .4e % .4e ",vx_t[current_index], vy_t[current_index], vz_t[current_index]);
1580 if (outputflags[4])
1581 fprintf(ptFile,"% .4e ",temp_t_m[current_index]);
1582 if (outputflags[5])
1583 fprintf(ptFile,"% .4e ",age_t[current_index]);
1584 fprintf(ptFile,"\n");
1585 }
1586 if (BLANK) fprintf(ptFile,"\n\n");
1587 }
1588 }
1589 }
1590 fprintf(ptFile,"# --------- cold layer -----------------------------------------------------------------------\n");
1591 for (k=0;k<kcmax+1;++k){
1592 for (j=0;j<jmax+1;++j){
1593 for (i=0;i<imax+1;++i){
1594 if( (outputflags[1]) && !((maske[j*(imax+1) + i]==0) || (maske[j*(imax+1) + i]==3)) ) continue;
1595 current_index = k*elements_in_layer + j*(imax+1) + i;
1596 if (outputflags[2])
1597 fprintf(ptFile,"%+.4e %+.4e %+.4e ",xi[i], eta[j], z_c[current_index]);
1598 if (outputflags[3])
1599 fprintf(ptFile,"% .4e % .4e % .4e ",vx_c[current_index], vy_c[current_index], vz_c[current_index]);
1600 if (outputflags[4])
1601 fprintf(ptFile,"% .4e ",temp_c[current_index]);
1602 if (outputflags[5])
1603 fprintf(ptFile,"% .4e ",age_c[current_index]);
1604 fprintf(ptFile,"\n");
1605 }
1606 if (BLANK) fprintf(ptFile,"\n\n\n");
1607 else fprintf(ptFile,"\n");
1608 }
1609 }
1610 }
1611 fclose(ptFile);
1612 }
1613 /* open file for writing of 2d variables */
1614 sprintf(filename,"%5s%2s_2d%s", runname, ergnum, suffix);
1615 if((ptFile=fopen(filename, "w"))==NULL){
1616 printf("\a Could not open file for writing of 3d ASCII-data\n");
1617 return;
1618 }else{
1619 printf("| writing on 2d output file %s\n", filename);
1620 fprintf(ptFile,"# *********************************************************\n");
1621 fprintf(ptFile,"# ASCCI output file of run %s, file no. %s\n", runname, ergnum);
1622 if (time(&how_late_is_it)==-1){
1623 printf("WARNING: Could not evaluate current time\n");
1624 fprintf(ptFile,"# no specific date was able to be inquiered\n");
1625 }else{
1626 fprintf(ptFile,"# written on %s", ctime(&how_late_is_it));
1627 }
1628
1629 #if defined(__APPLE__) || defined(MINGW32) || defined(WIN32) || defined(BSD)
1630 if ( 1 )
1631 #else
1632 if (cuserid(user_name)==NULL)
1633 #endif
1634 fprintf(ptFile,"# no user id was able to be inquiered\n");
1635 else
1636 fprintf(ptFile,"# by %10s\n", user_name);
1637 fprintf(ptFile,"# *********************************************************\n");
1638 fprintf(ptFile,"# 2d variables:\n");
1639 fprintf(ptFile,"# ");
1640 if (outputflags[11])
1641 fprintf(ptFile,"x y ");
1642 if (outputflags[6])
1643 fprintf(ptFile,"b s ");
1644 if (outputflags[7])
1645 fprintf(ptFile,"Drainage ");
1646 if (outputflags[8])
1647 fprintf(ptFile,"glaciation ");
1648 if (outputflags[9])
1649 fprintf(ptFile,"type_base ");
1650 if (outputflags[10])
1651 fprintf(ptFile,"flux_x flux_y ");
1652 fprintf(ptFile,"\n# --------------------------------------------------------------------------------------------------------\n");
1653 for (j=0;j<jmax+1;++j){
1654 for (i=0;i<imax+1;++i){
1655 current_index = j*(imax+1) + i;
1656 /* printf("%d ++ %d == %d, %d\n", outputflags[1], !((maske[current_index]==0) || (maske[current_index]==3)), ( (outputflags[1]==0) && !((maske[current_index]==0) || (maske[current_index]==3)) ), maske[current_index]); */
1657 if( (outputflags[1]==0) && !((maske[current_index]==0) || (maske[current_index]==3)) ) continue;
1658 if (outputflags[11])
1659 fprintf(ptFile,"% .4e % .4e", xi[i], eta[j]);
1660 if (outputflags[6]){
1661 /* printf("%d %d\n", i,j); */
1662 if (outputflags[0]){
1663 fprintf(ptFile,"% .4e % .4e", z_t[current_index], z_c[kcmax*elements_in_layer + current_index]);
1664 }
1665 else
1666 fprintf(ptFile,"% .4e % .4e", z_c[current_index], z_c[kcmax*elements_in_layer + current_index]);
1667 }
1668 if (outputflags[7])
1669 fprintf(ptFile,"% .4e ",Q_tld[current_index]);
1670 if (outputflags[8])
1671 fprintf(ptFile,"% 1d ",maske[current_index]);
1672 if (outputflags[9])
1673 fprintf(ptFile,"% 1d ",n_cts[current_index]);
1674 if (outputflags[10])
1675 fprintf(ptFile,"% .4e % .4e",qx[current_index],qy[current_index]);
1676 fprintf(ptFile,"\n");
1677 }
1678 if (BLANK) fprintf(ptFile,"\n\n");
1679 }
1680 fclose(ptFile);
1681 }
1682 printf("---------------------------------------------------------------\n");
1683 }
1684 /*******************************************************
1685 calculate staggered grid
1686 ********************************************************/
get_staggered_grid(float * xi,float * eta,float * z_in,int imax,int jmax,int kmax,float * deltaX,float * staggered_grid)1687 int get_staggered_grid(float *xi, /* unscaled x coordinate index i: 0,imax */
1688 float *eta, /* unscaled y coordinate index j from 0,jmax */
1689 float *z_in, /* unscaled z coordinate index i: 0,imax, j: 0,jmax, kc: 0,kcmax */
1690 int imax, /* grid steps in xi-direction */
1691 int jmax, /* grid steps in eta-direction */
1692 int kmax, /* grid steps in z-direction in cold ice layer */
1693 float *deltaX, /* stepsize of original grid */
1694 float *staggered_grid){ /* coords of staggered FEM grid */
1695
1696 register int i,j,k,n;
1697 int nodes_in_layer_of_staggered_grid, nodes_in_staggered_grid, auxiliary;
1698 float delta_x, *x,*y,*z;
1699
1700 nodes_in_layer_of_staggered_grid = (imax+2)*(jmax+2);
1701 nodes_in_staggered_grid =nodes_in_layer_of_staggered_grid*(kmax+1);
1702
1703 x=(float *) malloc((size_t) (imax+2)*sizeof(float));
1704 y=(float *) malloc((size_t) (jmax+2)*sizeof(float));
1705 z=(float *) malloc((size_t) nodes_in_staggered_grid*sizeof(float));
1706
1707 delta_x= 0.5*(*deltaX);
1708
1709 for (i=0;i<imax+1;++i)
1710 x[i]=xi[i] - delta_x;
1711 x[imax+1]=xi[imax] + delta_x;
1712 for (j=0;j<jmax+1;++j)
1713 y[j]=eta[j] - delta_x;
1714 y[jmax+1]=eta[jmax] + delta_x;
1715
1716 auxiliary = get_interpolated_property_on_staggered_grid(imax, jmax, kmax, z_in, z);
1717 if (auxiliary != nodes_in_staggered_grid){
1718 printf("WARNING: Returned value of interpolated z-values %d does not match calculated %d\n", auxiliary, nodes_in_staggered_grid);
1719 free(x);free(y);free(z);
1720 return -1;
1721 }
1722 for (k=0,n=0;k<kmax+1;++k){
1723 for(j=0;j<jmax+2;++j){
1724 for (i=0;i<imax+2;++i,++n){
1725 staggered_grid[(k*nodes_in_layer_of_staggered_grid + j*(imax+2) + i)*3]=x[i];
1726 staggered_grid[(k*nodes_in_layer_of_staggered_grid + j*(imax+2) + i)*3+1]=y[j];
1727 staggered_grid[(k*nodes_in_layer_of_staggered_grid + j*(imax+2) + i)*3+2]=z[k*nodes_in_layer_of_staggered_grid + j*(imax+2) + i];
1728 /* printf("%f %f %f\n", staggered_grid[(k*nodes_in_layer_of_staggered_grid + j*(imax+2) + i)*3], staggered_grid[(k*nodes_in_layer_of_staggered_grid + j*(imax+2) + i)*3+1], staggered_grid[(k*nodes_in_layer_of_staggered_grid + j*(imax+2) + i)*3 + 2]); */
1729 }
1730 }
1731 /* printf("\n"); */
1732 }
1733 free(x);free(y);free(z);
1734 return n;
1735 }
1736
1737 /*******************************************************
1738 interpolate property on staggered grid
1739 ********************************************************/
get_interpolated_property_on_staggered_grid(int imax,int jmax,int kmax,float * property_in,float * property_out)1740 int get_interpolated_property_on_staggered_grid(int imax,
1741 int jmax,
1742 int kmax,
1743 float *property_in,
1744 float *property_out){
1745 register int i,j,k,n;
1746 int nodes_in_layer_of_staggered_grid, nodes_in_layer;
1747
1748 nodes_in_layer = (imax+1)*(jmax+1);
1749 nodes_in_layer_of_staggered_grid = (imax+2)*(jmax+2);
1750
1751 /* inside array */
1752 for (k=0,n=0;k<kmax+1;++k){
1753 for (j=1;j<jmax+1;++j){
1754 for (i=1;i<imax+1;++i){
1755 property_out[k*nodes_in_layer_of_staggered_grid + j*(imax+2) + i]=0.25*(property_in[k*nodes_in_layer + j*(imax+1) + i-1] + property_in[k*nodes_in_layer + (j-1)*(imax+1) + i-1] + property_in[k*nodes_in_layer + j*(imax+1) + i] + property_in[k*nodes_in_layer + (j-1)*(imax+1) + i]);
1756 ++n;
1757 }
1758 }
1759 }
1760 /* at borders */
1761 for (k=0;k<kmax+1;++k){
1762 for (j=1;j<jmax+1;++j){
1763 property_out[k*nodes_in_layer_of_staggered_grid + j*(imax+2)]=0.5*(property_in[k*nodes_in_layer + j*(imax+1)] + property_in[k*nodes_in_layer + (j-1)*(imax+1)]);
1764 ++n;
1765 property_out[k*nodes_in_layer_of_staggered_grid + j*(imax+2) + imax+1]=0.5*(property_in[k*nodes_in_layer + j*(imax+1) + imax] + property_in[k*nodes_in_layer + (j-1)*(imax+1) + imax]);
1766 ++n;
1767 }
1768 for (i=1;i<imax+1;++i){
1769 property_out[k*nodes_in_layer_of_staggered_grid + i]= 0.5*(property_in[k*nodes_in_layer + i] + property_in[k*nodes_in_layer + i-1]);
1770 ++n;
1771 property_out[k*nodes_in_layer_of_staggered_grid + (jmax+1)*(imax+2) + i]=0.5*(property_in[k*nodes_in_layer + jmax*(imax+1) + i] + property_in[k*nodes_in_layer + jmax*(imax+1) + i-1]);
1772 ++n;
1773 }
1774 }
1775 /* at corners */
1776 for (k=0;k<kmax+1;++k){
1777 property_out[k*nodes_in_layer_of_staggered_grid] = property_in[k*nodes_in_layer];
1778 ++n;
1779 property_out[k*nodes_in_layer_of_staggered_grid + imax+1] = property_in[k*nodes_in_layer+ imax];
1780 ++n;
1781 property_out[k*nodes_in_layer_of_staggered_grid + (jmax+1)*(imax+2)] = property_in[k*nodes_in_layer+jmax*(imax+1)];
1782 ++n;
1783 property_out[k*nodes_in_layer_of_staggered_grid + j*(imax+2) + i] = property_in[k*nodes_in_layer+jmax*(imax+1)+imax];
1784 ++n;
1785 }
1786 return n;
1787 }
1788 /************************************************
1789 get information on glaciation of ice-sheet(s)
1790 *************************************************/
get_glaciation_info(int imax,int jmax,int * iced,int * maske)1791 int get_glaciation_info(int imax,
1792 int jmax,
1793 int *iced,
1794 int *maske){
1795 register int i,j,n;
1796
1797 for (j=0, n=0;j<jmax+1;++j){
1798 for (i=0;i<imax+1;++i){
1799 if ((maske[j*(imax+1) + i]==0) || (maske[j*(imax+1) + i]==3)){
1800 iced[j*(imax+1) + i]=n;
1801 ++n;
1802 }else{
1803 iced[j*(imax+1) + i]=-1;
1804 }
1805 }
1806 }
1807 return n;
1808 }
1809 /************************************************
1810 get 2d glaciation map for staggered grid
1811 *************************************************/
get_grid_map(int imax,int jmax,int * iced,int * gridmap)1812 int get_grid_map(int imax,
1813 int jmax,
1814 int *iced,
1815 int *gridmap){
1816 register int i,j,n;
1817 int nodes_in_layer_of_staggered_grid, nodes_in_layer, *staggered_iced;
1818 nodes_in_layer = (imax+1)*(jmax+1);
1819 nodes_in_layer_of_staggered_grid = (imax+2)*(jmax+2);
1820 if ((staggered_iced = (int *) malloc((size_t) nodes_in_layer_of_staggered_grid*sizeof(int)))==NULL){
1821 printf("ERROR during allocation of memory for glaciation information of staggered grid\n");
1822 return -1;
1823 }
1824 for (j=0;j<jmax+2;++j){
1825 for (i=0;i<imax+2;++i){
1826 staggered_iced[j*(imax+2) + i]=0;
1827 }
1828 }
1829 for (j=0;j<jmax+1;++j){
1830 for (i=0;i<imax+1;++i){
1831 if (iced[j*(imax+1) + i]!=-1){
1832 staggered_iced[j*(imax+2) + i]=1;
1833 staggered_iced[(j+1)*(imax+2) + i]=1;
1834 staggered_iced[j*(imax+2) + i+1]=1;
1835 staggered_iced[(j+1)*(imax+2) + i+1]=1;
1836 }
1837 }
1838 }
1839 for (j=0, n=0;j<jmax+2;++j){
1840 for (i=0;i<imax+2;++i){
1841 if (staggered_iced[j*(imax+2) + i]){
1842 gridmap[j*(imax+2) + i]=n;
1843 ++n;
1844 }else{
1845 gridmap[j*(imax+2) + i]=-1;
1846 }
1847 }
1848 }
1849 free(staggered_iced);
1850 return n;
1851 }
1852 /**************************************************
1853 get 2d information on boundaries of ice-sheet(s)
1854 ***************************************************/
get_glaciation_boundary_info(int imax,int jmax,int * iced,int * boundary)1855 int get_glaciation_boundary_info(int imax,
1856 int jmax,
1857 int *iced,
1858 int *boundary){
1859 register int i,j,n,nn;
1860 /* inside array */
1861 for (j=1,n=0,nn=0;j<jmax;++j){
1862 for (i=1;i<imax;++i){
1863 if (iced[j*(imax+1) + i]!=-1){
1864 if (iced[j*(imax+1) + i+1]==-1){ /* east */
1865 boundary[iced[j*(imax+1) + i]*4+0]=1;
1866 ++n;
1867 }else boundary[iced[j*(imax+1) + i]*4+0]=0;
1868 if (iced[(j+1)*(imax+1) + i]==-1){ /* north */
1869 boundary[iced[j*(imax+1) + i]*4+1]=1;
1870 ++n;
1871 }else boundary[iced[j*(imax+1) + i]*4+1]=0;
1872 if (iced[j*(imax+1) + i-1]==-1){ /* west */
1873 boundary[iced[j*(imax+1) + i]*4+2]=1;
1874 ++n;
1875 }else boundary[iced[j*(imax+1) + i]*4+2]=0;
1876 if (iced[(j-1)*(imax+1) + i]==-1){ /* south */
1877 boundary[iced[j*(imax+1) + i]*4+3]=1;
1878 ++n;
1879 }else boundary[iced[j*(imax+1) + i]*4+3]=0;
1880 }
1881 }
1882 }
1883 /* at borders */
1884 for (j=1;j<jmax;++j){
1885 /* most western row (i==0)*/
1886 if (iced[j*(imax+1)]!=-1){
1887 boundary[iced[j*(imax+1)]*4+2]=1; /* west is always boundary */
1888 ++n;
1889 if (iced[j*(imax+1) +1]==-1){ /* east */
1890 boundary[iced[j*(imax+1)]*4+0]=1;
1891 ++n;
1892 }else boundary[iced[j*(imax+1)]*4+0]=0;
1893 if (iced[(j+1)*(imax+1)]==-1){ /* north */
1894 boundary[iced[j*(imax+1)]*4+1]=1;
1895 ++n;
1896 }else boundary[iced[j*(imax+1)]*4+1]=0;
1897 if (iced[(j-1)*(imax+1)]==-1){ /* south */
1898 boundary[iced[j*(imax+1)]*4+3]=1;
1899 ++n;
1900 }else boundary[iced[j*(imax+1)]*4+3]=0;
1901 }
1902 /* most eastern row (i==imax)*/
1903 if (iced[j*(imax+1)+imax]!=-1){
1904 boundary[iced[j*(imax+1) + imax]*4+0]=1;/* east is always boundary */
1905 ++n;
1906 if (iced[(j+1)*(imax+1) + imax]==-1){ /* north */
1907 boundary[iced[j*(imax+1) + imax]*4+1]=1;
1908 ++n;
1909 }else boundary[iced[j*(imax+1) + imax]*4+1]=0;
1910 if (iced[j*(imax+1) + imax-1]==-1){ /* west */
1911 boundary[iced[j*(imax+1) + imax]*4+2]=1;
1912 ++n;
1913 }else boundary[iced[j*(imax+1) + imax]*4+2]=0;
1914 if (iced[(j-1)*(imax+1) + imax]==-1){ /* south */
1915 boundary[iced[j*(imax+1) + imax]*4+3]=1;
1916 ++n;
1917 }else boundary[iced[j*(imax+1) + imax]*4+3]=0;
1918 }
1919 }
1920 for (i=1;i<imax;++i){
1921 /* most northern row (j==jmax)*/
1922 if (iced[jmax*(imax+1) + i]!=-1){
1923 boundary[iced[jmax*(imax+1) + i]*4+1]=1;/* north is always boundary */
1924 ++n;
1925 if (iced[jmax*(imax+1) + i+1]==-1){ /* east */
1926 boundary[iced[jmax*(imax+1) + i]*4+0]=1;
1927 ++n;
1928 }else boundary[iced[jmax*(imax+1) + i]*4+0]=0;
1929 if (iced[jmax*(imax+1) + i-1]==-1){ /* west */
1930 boundary[iced[jmax*(imax+1) + i]*4+2]=1;
1931 ++n;
1932 }else boundary[iced[jmax*(imax+1) + i]*4+2]=0;
1933 if (iced[(jmax-1)*(imax+1) + i]==-1){ /* south */
1934 boundary[iced[jmax*(imax+1) + i]*4+3]=1;
1935 ++n;
1936 }else boundary[iced[jmax*(imax+1) + i]*4+3]=0;
1937 }
1938 /* most southern (j==0)*/
1939 if (iced[i]!=-1){
1940 boundary[iced[i]*4+3]=1; /* south is always boundary */
1941 ++n;
1942 if (iced[i+1]==-1){ /* east */
1943 boundary[iced[i]*4+0]=1;
1944 ++n;
1945 }else boundary[iced[i]*4+0]=0;
1946 if (iced[1*(imax+1) + i]==-1){ /* north */
1947 boundary[iced[i]*4+1]=1;
1948 ++n;
1949 }else boundary[iced[i]*4+1]=0;
1950 if (iced[i-1]==-1){ /* west */
1951 boundary[iced[i]*4+2]=1;
1952 ++n;
1953 }else boundary[iced[i]*4+2]=0;
1954 }
1955 }
1956 /* at corners */
1957 /* northeast */
1958 if (iced[jmax*(imax+1) + imax]!=-1){
1959 boundary[iced[jmax*(imax+1) + imax]*4+0]=1; /* east is always boundary */
1960 ++n;
1961 boundary[iced[jmax*(imax+1) + i]*4+1]=1;/* north is always boundary */
1962 ++n;
1963 if (iced[jmax*(imax+1) + imax-1]==-1){ /* west */
1964 boundary[iced[jmax*(imax+1) + imax]*4+2]=1;
1965 ++n;
1966 }else boundary[iced[jmax*(imax+1) + imax]*4+2]=0;
1967 if (iced[(jmax-1)*(imax+1) + imax]==-1){ /* south */
1968 boundary[iced[jmax*(imax+1) + imax]*4+3]=1;
1969 ++n;
1970 }else boundary[iced[jmax*(imax+1) + imax]*4+3]=0;
1971 }
1972 /* northwest */
1973 if (iced[jmax*(imax+1)]!=-1){
1974 if (iced[jmax*(imax+1)+1]==-1){ /* east */
1975 boundary[iced[jmax*(imax+1)]*4+0]=1;
1976 ++n;
1977 }else boundary[iced[jmax*(imax+1)]*4+0]=0;
1978 boundary[iced[jmax*(imax+1)]*4+1]=1;/* north is always boundary */
1979 ++n;
1980 boundary[iced[jmax*(imax+1)]*4+2]=1;/* west is always boundary */
1981 ++n;
1982 if (iced[(jmax-1)*(imax+1)]==-1){ /* south */
1983 boundary[iced[jmax*(imax+1)]*4+3]=1;
1984 ++n;
1985 }else boundary[iced[jmax*(imax+1)]*4+3]=0;
1986 }
1987 /* southwest */
1988 if (iced[0]!=-1){
1989 if (iced[1]==-1){ /* east */
1990 boundary[iced[0]*4+0]=1;
1991 ++n;
1992 }else boundary[iced[0]*4+0]=0;
1993 if (iced[1*(imax+1)]==-1){ /* north */
1994 boundary[iced[0]*4+1]=1;
1995 ++n;
1996 }else boundary[iced[0]*4+1]=0;
1997 boundary[iced[0]*4+2]=1;/* west is always boundary */
1998 ++n;
1999 boundary[iced[0]*4+3]=1;/* south is always boundary */
2000 ++n;
2001 }
2002 /* southeast */
2003 if (iced[imax]!=-1){
2004 boundary[iced[imax]*4+0]=1;/* east is always boundary */
2005 ++n;
2006 if (iced[1*(imax+1) + imax]==-1){ /* north */
2007 boundary[iced[imax]*4+1]=1;
2008 ++n;
2009 }else boundary[iced[imax]*4+1]=0;
2010 if (iced[imax-1]==-1){ /* west */
2011 boundary[iced[imax]*4+2]=1;
2012 ++n;
2013 }else boundary[iced[imax]*4+2]=0;
2014 boundary[iced[imax]*4+3]=1;/* south is always boundary */
2015 ++n;
2016 }
2017 return n;
2018 }
2019 /*******************************************************
2020 transform integer into float array
2021 ********************************************************/
make_float_from_integer_scalar_field(int * input_property,float * output_property,int number_of_nodes,int reorder_ice_land_sea_mask)2022 void make_float_from_integer_scalar_field(int *input_property,
2023 float *output_property,
2024 int number_of_nodes,
2025 int reorder_ice_land_sea_mask){
2026 register int n;
2027 if (reorder_ice_land_sea_mask){
2028 for (n=0;n<number_of_nodes;++n){
2029 if ((input_property[n]==0) || (input_property[n]==3)) output_property[n] = 0; /* ice */
2030 else if (input_property[n]==1) output_property[n] = 1; /* land */
2031 else output_property[n] = 2; /* sea */
2032 }
2033 }else{
2034 for (n=0;n<number_of_nodes;++n)
2035 output_property[n] = (float) input_property[n];
2036 }
2037 return;
2038 }
2039
2040
2041
2042
2043
2044
2045 /****************************************************************
2046 Output of SICOPOLIS grid (not staggerd) in ELMER Solver format
2047 *****************************************************************/
FC_FUNC(solvergrid,SOLVERGRID)2048 void STDCALLBULL FC_FUNC(solvergrid,SOLVERGRID)(float *xi, /* unscaled x coordinate index i: 0,imax */
2049 float *eta, /* unscaled y coordinate index j from 0,jmax */
2050 float *z_c, /* unscaled z coordinate index i: 0,imax, j: 0,jmax, kc: 0,kcmax */
2051 float *z_t, /* unscaled y coordinate index i: 0,imax, j: 0,jmax, kt: 0,kcmax */
2052 int *imax_in, /* grid steps in xi-direction */
2053 int *jmax_in, /* grid steps in eta-direction */
2054 int *kcmax_in, /* grid steps in z-direction in cold ice layer */
2055 int *ktmax_in, /* grid steps in z-direction in temperate ice layer */
2056 FC_CHAR_PTR(runname,runname_l), /*name of run*/
2057 FC_CHAR_PTR(ergnum,ergnum_l), /*number of file*/
2058 int *maske, /*mask of vertex type */
2059 float *deltaX, /* stepsize of grid */
2060 int *flag)
2061 {
2062 register int i, j, k, l, m, n;
2063 int number_of_bulk_elements, number_of_boundary_elements, number_of_iced_collums, elements_in_one_layer, number_of_nodes, nodes_of_element[8], *nodes_of_side_element, *parent_of_side_element, nodes_in_one_layer, *iced, idnr, parent_element, sidebulk=0;
2064 int imax, jmax, kcmax, ktmax, kmax = 0;
2065 float min_max_xyz[2][3], *bottom, freesurf, *delta_z;
2066 float actual_scaled_coord[3];
2067 char groupid[4], filename[80], yes_no;
2068 FILE *ptFile;
2069
2070 /* constants */
2071 imax= *imax_in;
2072 jmax= *jmax_in;
2073 kcmax= *kcmax_in;
2074 ktmax= *ktmax_in;
2075
2076 /* print out little summary */
2077 printf("---------------------------------------------------------------\n");
2078 printf("| Output of SICOPOLIS Grid for ELMER Solver\n");
2079 printf("---------------------------------------------------------------\n");
2080 printf("| imax/jmax/kcmax/ktmax=%d/%d/%d/%d\n",imax, jmax, kcmax, ktmax);
2081 printf("---------------------------------------------------------------\n");
2082 printf("| nodes in original grid:\n");
2083 printf("| cold layer: %d=%d * %d\n", (imax+1)*(jmax+1)*(kcmax+1), (imax+1)*(jmax+1), (kcmax+1));
2084 printf("| temperate layer: %d=%d * (%d - 1)\n", (imax+1)*(jmax+1)*(ktmax), (imax+1)*(jmax+1), (ktmax+1));
2085 printf("| -------------\n");
2086 printf("| %d=%d * %d\n", (imax+1)*(jmax+1)*(ktmax+1+kcmax), (imax+1)*(jmax+1), ktmax+1+kcmax);
2087 printf("---------------------------------------------------------------\n");
2088 printf("| elements in original grid:\n");
2089 printf("| cold layer: %d\n", (imax)*(jmax)*kcmax);
2090 printf("| temperate layer: %d\n", (*flag)*(imax)*(jmax)*(ktmax));
2091 printf("---------------------------------------------------------------\n");
2092 printf("| x = %.1f -> %.1f , y = %.1f -> %.1f, dx=%.1f\n", xi[0], xi[imax], eta[0], eta[jmax], *deltaX);
2093 printf("---------------------------------------------------------------\n");
2094
2095 /* inquire number of grid levels */
2096 /* ------------------------------*/
2097 while (kmax < 1){
2098 printf("| How many vertical grid layers the grid shall contain? \n");
2099 scanf("%d", &kmax);
2100 printf("| \n");
2101 if (kmax < 1) printf("| No. of vertical (current value %d) levels must exceed 1!\n", kmax);
2102 }
2103 printf("Thanks. Taking %d vertical layers\n", kmax);
2104 printf("---------------------------------------------------------------\n");
2105
2106 /* inquire if sidebulk on free surface is desired */
2107 /* -----------------------------------------------*/
2108 printf("| sidebulk elements on free surface (y/n)? \n");
2109 scanf("%1s", &yes_no);
2110 printf("| \n");
2111 if ((yes_no == 'y') || (yes_no == 'Y') ){
2112 printf("| Thanks. Writing sidebulk elements on free surface\n");
2113 sidebulk = 1;
2114 }else{
2115 printf("| Thanks. No output of sidebulk elements on free surface\n");
2116 sidebulk = 0;
2117 }
2118 printf("---------------------------------------------------------------\n");
2119
2120 /* calculate constants for output mesh*/
2121 /* -----------------------------------*/
2122 nodes_in_one_layer = (imax+1)*(jmax+1);
2123 elements_in_one_layer = (imax)*(jmax);
2124 number_of_nodes = nodes_in_one_layer*(kmax+1);
2125 number_of_bulk_elements = elements_in_one_layer*kmax;
2126 number_of_boundary_elements = 2*elements_in_one_layer + 2*(imax+jmax)*kmax;
2127
2128 /* allocate arrays */
2129 /* ----------------*/
2130 iced = (int *) malloc((size_t) (imax+1)*(jmax+1)*sizeof(int));
2131 if (iced == NULL){
2132 printf("ERROR in allocating memory for glaciation information\n");
2133 free(iced);
2134 return;
2135 }
2136 delta_z = (float *) malloc((size_t) (imax+1)*(jmax+1)*sizeof(float));
2137 if (delta_z == NULL){
2138 printf("ERROR in allocating memory for vertical grid step sizes\n");
2139 free(iced);free(delta_z);
2140 return;
2141 }
2142 bottom = (float *) malloc((size_t) (imax+1)*(jmax+1)*sizeof(float));
2143 if (bottom == NULL){
2144 printf("ERROR in allocating memory for bottom topography\n");
2145 free(iced);free(delta_z);free(bottom);
2146 return;
2147 }
2148
2149 /* get glaciation info */
2150 /* --------------------*/
2151 number_of_iced_collums = get_glaciation_info(imax,jmax,iced,maske);
2152 if (number_of_iced_collums<1){
2153 printf("WARNING: domain seems to be ice-free\n");
2154 }
2155 printf("| number of iced colums: %d out of %d (%3.2f %%)\n", number_of_iced_collums, (imax+1)*(jmax+1), ((float) number_of_iced_collums)*100.0/((float) (imax+1)*(jmax+1)));
2156 printf("---------------------------------------------------------------\n");
2157 printf("| number of nodes: %d=%d * %d\n",number_of_nodes, nodes_in_one_layer, kmax+1);
2158 printf("| number of bulk elements: %d=%d * %d\n",number_of_bulk_elements, elements_in_one_layer, kmax);
2159 printf("| number of sidebulk elements: %d\n", sidebulk*elements_in_one_layer);
2160 printf("| number of boundary elements: %d=2*[%d + (%d + %d) *%d] + %d\n", number_of_boundary_elements, elements_in_one_layer, imax, jmax, kmax, 2*(imax+jmax)*sidebulk);
2161 printf("| sidebulk boundary elements: %d=%d*2*(%d + %d)\n", sidebulk*2*(imax+jmax), sidebulk, imax, jmax);
2162 printf("---------------------------------------------------------------\n");
2163
2164 /* writing header file */
2165 /* --------------------*/
2166 sprintf(filename,"mesh.header");
2167 printf("| Writing mesh header file %s.\n",filename);
2168 if((ptFile=fopen(filename, "w"))==NULL){
2169 printf("\a Could not open Elmer mesh-file file for writing!\n");
2170 free(iced);free(delta_z);free(bottom);
2171 return;
2172 }
2173 fprintf(ptFile, "%d %d %d\n", number_of_nodes, number_of_bulk_elements + sidebulk*elements_in_one_layer, number_of_boundary_elements + sidebulk*2*(imax+jmax));
2174 fprintf(ptFile,"%d\n",2+sidebulk);
2175 if (sidebulk)
2176 fprintf(ptFile,"202 %d\n", 2*(imax+jmax));
2177 fprintf(ptFile,"404 %d\n", number_of_boundary_elements + sidebulk*elements_in_one_layer);
2178 fprintf(ptFile,"808 %d\n", number_of_bulk_elements);
2179
2180 printf("| succeeded in writting mesh header file %s.\n",filename);
2181 printf("---------------------------------------------------------------\n");
2182 fclose(ptFile);
2183
2184 /* init min/max */
2185 /* -------------*/
2186 min_max_xyz[0][2]= min_max_xyz[1][2] = z_c[0];
2187 min_max_xyz[0][0]= xi[0];
2188 min_max_xyz[1][0] = xi[imax];
2189 min_max_xyz[0][1]= eta[0];
2190 min_max_xyz[1][1] = eta[jmax];
2191
2192 /* get delta_z for specific i,j value */
2193 /* -----------------------------------*/
2194 for (j=0;j<jmax+1;++j){ /* loop over all j-values */
2195 for (i=0;i<imax+1;++i) { /* loop over all i-values */
2196 freesurf = z_c[kcmax*nodes_in_one_layer + j*(imax+1) + i];
2197 if (*flag){
2198 bottom[j*(imax+1) + i] = z_t[j*(imax+1) + i];
2199 }else{
2200 bottom[j*(imax+1) + i] = z_c[j*(imax+1) + i];
2201 }
2202 if (min_max_xyz[0][2]>bottom[j*(imax+1) + i]) min_max_xyz[0][2]=bottom[j*(imax+1) + i];
2203 if (min_max_xyz[1][2]<freesurf) min_max_xyz[0][2]=freesurf;
2204 delta_z[j*(imax+1) + i] = (freesurf - bottom[j*(imax+1) + i])/((float) kmax);
2205 /* printf("dz=%f = %f - %f \n", delta_z[j*(imax+1) + i], freesurf, bottom[j*(imax+1) + i]); */
2206 /* printf("%d %d % f\n", i,j,delta_z[j*(imax+1) + i]); */
2207 if (delta_z[j*(imax+1) + i] < 0.0){
2208 printf("\a delta z = %f - %f for (%d, %d) = %f < 0\n", freesurf, bottom[j*(imax+1) + i], i, j, delta_z[j*(imax+1) + i]);
2209 free(iced);free(delta_z);free(bottom);
2210 fclose(ptFile);
2211 return;
2212 }else if((delta_z[j*(imax+1) + i] < MINHEIGHT/((float) kmax))|| ((maske[j*(imax+1) + i]!=0) && (maske[j*(imax+1) + i]!=3))){
2213 delta_z[j*(imax+1) + i] = MINHEIGHT/((float) kmax);
2214 }
2215 }
2216 }
2217 printf("| succeeded in calculating grid step sizes %s.\n",filename);
2218 printf("---------------------------------------------------------------\n");
2219
2220 /* writing nodes file */
2221 /* -------------------*/
2222 sprintf(filename,"mesh.nodes");
2223 printf("| Writing node-file %s.\n",filename);
2224 if((ptFile=fopen(filename, "w"))==NULL){
2225 printf("\a Could not open Elmer mesh-file file for writing!\n");
2226 free(iced);free(delta_z);free(bottom);
2227 return;
2228 }
2229 for (n=0,k=0;k<kmax+1;++k){ /* loop over all levels */
2230 for (j=0;j<jmax+1;++j){ /* loop over all j-values */
2231 for (i=0;i<imax+1;++i) { /* loop over all i-values */
2232 ++n;
2233 fprintf(ptFile, "%d -1 %.8f %.8f %.8f\n", n, eta[i], xi[j], bottom[j*(imax+1) + i] + delta_z[j*(imax+1) + i] * ((float) k));
2234 }
2235 }
2236 }
2237 fclose(ptFile);
2238 if (n != number_of_nodes){
2239 printf("\a %d written nodes does not match %d calculated nodes in grid\n", n, number_of_nodes);
2240 return;
2241 }else{
2242 printf("| succeeded in writing node file %s.\n",filename);
2243 printf("---------------------------------------------------------------\n");
2244 }
2245
2246 /* writing element file */
2247 /* ---------------------*/
2248 sprintf(filename,"mesh.elements");
2249 printf("| Writing bulk element file %s.\n",filename);
2250 if((ptFile=fopen(filename, "w"))==NULL){
2251 printf("\a Could not open Elmer mesh-file file for writing!\n");
2252 free(iced);free(delta_z);free(bottom);
2253 return;
2254 }
2255 /* write bulk elements */
2256 idnr = 1;
2257 for (n=0,k=0;k<kmax;++k){ /* loop over all levels */
2258 for (j=0;j<jmax;++j){ /* loop over all j-values */
2259 for (i=0;i<imax;++i) { /* loop over all i-values */
2260 n++;
2261 for (l=0; l<2; ++l){/* lower level: n=0; upper level n=1; each counterclkws */
2262 nodes_of_element[l*4] = (k+l)*nodes_in_one_layer + j*(imax+1) + i;
2263 nodes_of_element[l*4+1] = (k+l)*nodes_in_one_layer + j*(imax+1) + i+1;
2264 nodes_of_element[l*4+2] = (k+l)*nodes_in_one_layer + (j+1)*(imax+1) + i+1;
2265 nodes_of_element[l*4+3] = (k+l)*nodes_in_one_layer + (j+1)*(imax+1) + i;
2266 }
2267 fprintf(ptFile,"%d %d 808 %d %d %d %d %d %d %d %d\n", n, idnr,
2268 nodes_of_element[0]+1, nodes_of_element[1]+1, nodes_of_element[2]+1, nodes_of_element[3]+1,
2269 nodes_of_element[4]+1, nodes_of_element[5]+1, nodes_of_element[6]+1, nodes_of_element[7]+1);
2270 }
2271 }
2272 }
2273 /* write sidebulk elements (if demanded) */
2274 if (sidebulk){
2275 idnr = 2;
2276 for (m=0,j=0;j<jmax;++j){ /* loop over all j-values */
2277 for (i=0;i<imax;++i) { /* loop over all i-values */
2278 ++m;
2279 nodes_of_element[0] = kmax*nodes_in_one_layer + j*(imax+1) + i;
2280 nodes_of_element[1] = kmax*nodes_in_one_layer + j*(imax+1) + i+1;
2281 nodes_of_element[2] = kmax*nodes_in_one_layer + (j+1)*(imax+1) + i+1;
2282 nodes_of_element[3] = kmax*nodes_in_one_layer + (j+1)*(imax+1) + i;
2283 fprintf(ptFile,"%d %d 404 %d %d %d %d\n", n+m, idnr,
2284 nodes_of_element[0]+1, nodes_of_element[1]+1, nodes_of_element[2]+1, nodes_of_element[3]+1);
2285 }
2286 }
2287 }else
2288 m=0;
2289 fclose(ptFile);
2290 if (n+m != number_of_bulk_elements + sidebulk*elements_in_one_layer){
2291 printf("\a %d=%d + %d written bulk elements does not match %d= %d + %d calculated elements in grid\n",
2292 n+m, n,m, number_of_bulk_elements+elements_in_one_layer, number_of_bulk_elements, elements_in_one_layer);
2293 return;
2294 }else{
2295 printf("| %d number of bulk elements written.\n", n);
2296 printf("| %d number of sidebulk elements written.\n", m);
2297 printf("| succeeded in writting bulk element file %s.\n",filename);
2298 printf("---------------------------------------------------------------\n");
2299 }
2300
2301 /* writing boundary element file */
2302 /* ------------------------------*/
2303 sprintf(filename,"mesh.boundary");
2304 printf("| Writing boundary element file %s.\n",filename);
2305 if((ptFile=fopen(filename, "w"))==NULL){
2306 printf("\a Could not open Elmer boundary element file for writing!\n");
2307 free(iced);free(delta_z);free(bottom);
2308 return;
2309 }
2310 /* base */
2311 idnr = 1;
2312 printf("| lower (base, z=0) boundary for bulk; idnr=%d\n", idnr);
2313 for (n=0,j=0;j<jmax;++j){ /* loop over all j-values */
2314 for (i=0;i<imax;++i) { /* loop over all i-values */
2315 n++;
2316 nodes_of_element[0] = j*(imax+1) + i;
2317 nodes_of_element[1] = j*(imax+1) + i+1;
2318 nodes_of_element[2] = (j+1)*(imax+1) + i+1;
2319 nodes_of_element[3] = (j+1)*(imax+1) + i;
2320 parent_element = j*imax + i;
2321 if (parent_element+1 > number_of_bulk_elements){
2322 printf("parent element %d > number of bulk elments %d\n",parent_element+1, number_of_bulk_elements);
2323 return;
2324 }
2325 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n",
2326 n, idnr, parent_element+1, nodes_of_element[0]+1, nodes_of_element[1]+1, nodes_of_element[2]+1, nodes_of_element[3]+1);
2327 }
2328 }
2329 /* free surface */
2330 idnr = 2;
2331 printf("| upper (free surface, z=max) boundary for bulk; idnr=%d\n", idnr);
2332 for (j=0;j<jmax;++j){ /* loop over all j-values */
2333 for (i=0;i<imax;++i) { /* loop over all i-values */
2334 n++;
2335 nodes_of_element[0] = kmax*nodes_in_one_layer + j*(imax+1) + i;
2336 nodes_of_element[1] = kmax*nodes_in_one_layer + j*(imax+1) + i+1;
2337 nodes_of_element[2] = kmax*nodes_in_one_layer + (j+1)*(imax+1) + i+1;
2338 nodes_of_element[3] = kmax*nodes_in_one_layer + (j+1)*(imax+1) + i;
2339 parent_element = (kmax-1)*elements_in_one_layer + j*imax + i; /* same numbering as in bulk */
2340 if (parent_element > number_of_bulk_elements){
2341 printf("parent element %d > number of bulk elments %d\n",parent_element+1, number_of_bulk_elements);
2342 return;
2343 }
2344 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n",
2345 n, idnr, parent_element+1, nodes_of_element[0]+1, nodes_of_element[1]+1, nodes_of_element[2]+1, nodes_of_element[3]+1);
2346 }
2347 }
2348 /* side faces: */
2349 /* south */
2350 idnr = 3;
2351 j=0;
2352 printf("| southern (y=0) boundary for bulk; idnr=%d\n", idnr);
2353 for (k=0;k<kmax;++k){ /* loop over all levels */
2354 for (i=0;i<imax;++i) { /* loop over all i-values */
2355 ++n;
2356 nodes_of_element[0] = k*nodes_in_one_layer + j*(imax+1) + i;
2357 nodes_of_element[1] = k*nodes_in_one_layer + j*(imax+1) + i+1;
2358 nodes_of_element[2] = (k+1)*nodes_in_one_layer + j*(imax+1) + i+1;
2359 nodes_of_element[3] = (k+1)*nodes_in_one_layer + j*(imax+1) + i;
2360 parent_element = k*elements_in_one_layer + j*imax + i; /* same numbering as in bulk */
2361 if (parent_element+1 > number_of_bulk_elements){
2362 printf("parent element %d > number of bulk elments %d\n",parent_element+1, number_of_bulk_elements);
2363 return;
2364 }
2365 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n",
2366 n, idnr, parent_element+1, nodes_of_element[0]+1, nodes_of_element[1]+1, nodes_of_element[2]+1, nodes_of_element[3]+1);
2367 }
2368 }
2369 /* west */
2370 idnr = 4;
2371 i=0;
2372 printf("| western (x=0) boundary for bulk; idnr=%d\n", idnr);
2373 for (k=0;k<kmax;++k){ /* loop over all levels */
2374 for (j=0;j<jmax;++j){ /* loop over all j-values */
2375 ++n;
2376 nodes_of_element[0] = k*nodes_in_one_layer + j*(imax+1) + i;
2377 nodes_of_element[1] = k*nodes_in_one_layer + (j+1)*(imax+1) + i;
2378 nodes_of_element[2] = (k+1)*nodes_in_one_layer + (j+1)*(imax+1) + i;
2379 nodes_of_element[3] = (k+1)*nodes_in_one_layer + j*(imax+1) + i;
2380 parent_element = k*elements_in_one_layer + j*imax + i; /* same numbering as in bulk */
2381 if (parent_element+1 > number_of_bulk_elements){
2382 printf("parent element %d > number of bulk elments %d\n",parent_element+1, number_of_bulk_elements);
2383 return;
2384 }
2385 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n",
2386 n, idnr, parent_element+1, nodes_of_element[0]+1, nodes_of_element[1]+1, nodes_of_element[2]+1, nodes_of_element[3]+1);
2387 }
2388 }
2389 /* north */
2390 idnr = 5;
2391 j=jmax;
2392 printf("| northern (y=max) boundary for bulk; idnr=%d\n", idnr);
2393 for (k=0;k<kmax;++k){ /* loop over all levels */
2394 for (i=0;i<imax;++i) { /* loop over all i-values */
2395 nodes_of_element[0] = k*nodes_in_one_layer + j*(imax+1) + i;
2396 nodes_of_element[1] = k*nodes_in_one_layer + j*(imax+1) + i+1;
2397 nodes_of_element[2] = (k+1)*nodes_in_one_layer + j*(imax+1) + i+1;
2398 nodes_of_element[3] = (k+1)*nodes_in_one_layer + j*(imax+1) + i;
2399 ++n;
2400 parent_element = k*elements_in_one_layer + (jmax-1)*imax + i; /* same numbering as in bulk */
2401 if (parent_element+1 > number_of_bulk_elements){
2402 printf("parent element %d > number of bulk elments %d\n",parent_element+1, number_of_bulk_elements);
2403 return;
2404 }
2405 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n",
2406 n, idnr, parent_element + 1, nodes_of_element[0] + 1, nodes_of_element[1] + 1, nodes_of_element[2] + 1, nodes_of_element[3] + 1);
2407 }
2408 }
2409 /* east */
2410 idnr = 6;
2411 i=imax;
2412 printf("| eastern (x=max) boundary for bulk; idnr=%d\n", idnr);
2413 for (k=0;k<kmax;++k){ /* loop over all levels */
2414 for (j=0;j<jmax;++j){ /* loop over all j-values */
2415 nodes_of_element[0] = k*nodes_in_one_layer + j*(imax+1) + i;
2416 nodes_of_element[1] = k*nodes_in_one_layer + (j+1)*(imax+1) + i;
2417 nodes_of_element[2] = (k+1)*nodes_in_one_layer + (j+1)*(imax+1) + i;
2418 nodes_of_element[3] = (k+1)*nodes_in_one_layer + j*(imax+1) + i;
2419 ++n;
2420 parent_element = k*elements_in_one_layer + j*imax + i-1; /* same numbering as in bulk */
2421 if (parent_element >= number_of_bulk_elements){
2422 printf("parent element %d > number of bulk elments %d\n",parent_element, number_of_bulk_elements);
2423 /* return; */
2424 }
2425 fprintf(ptFile,"%d %d %d 0 404 %d %d %d %d\n",
2426 n, idnr, parent_element + 1, nodes_of_element[0] + 1, nodes_of_element[1] + 1, nodes_of_element[2] + 1, nodes_of_element[3] + 1);
2427 }
2428 }
2429 m=0;
2430 if (sidebulk){
2431 /* frame of sidebulk (i.e. free surface): */
2432 /* south */
2433 ++idnr;
2434 j=0;
2435 printf("| southern (y=0) boundary for sidebulk; idnr=%d\n", idnr);
2436 for (i=0;i<imax;++i) { /* loop over all i-values */
2437 ++m;
2438 nodes_of_element[0] = kmax*nodes_in_one_layer + j*(imax+1) + i;
2439 nodes_of_element[1] = kmax*nodes_in_one_layer + j*(imax+1) + i+1;
2440 parent_element = number_of_bulk_elements + j*imax + i;
2441 fprintf(ptFile,"%d %d %d 0 202 %d %d\n",
2442 n+m, idnr, parent_element + 1, nodes_of_element[0] + 1, nodes_of_element[1] + 1);
2443 }
2444 /* west */
2445 ++idnr;
2446 i=0;
2447 printf("| western (x=0) boundary for sidebulk; idnr=%d\n", idnr);
2448 for (j=0;j<jmax;++j){ /* loop over all j-values */
2449 ++m;
2450 nodes_of_element[0] = kmax*nodes_in_one_layer + j*(imax+1) + i;
2451 nodes_of_element[1] = kmax*nodes_in_one_layer + (j+1)*(imax+1) + i;
2452 parent_element = number_of_bulk_elements + j*imax + i;
2453 fprintf(ptFile,"%d %d %d 0 202 %d %d\n",
2454 n+m, idnr, parent_element + 1, nodes_of_element[0] + 1, nodes_of_element[1] + 1);
2455 }
2456 /* north */
2457 ++idnr;
2458 j=jmax;
2459 printf("| northern (y=max) boundary for sidebulk; idnr=%d\n", idnr);
2460 for (i=0;i<imax;++i) { /* loop over all i-values */
2461 ++m;
2462 nodes_of_element[0] = (kmax-1)*nodes_in_one_layer + j*(imax+1) + i;
2463 nodes_of_element[1] = (kmax-1)*nodes_in_one_layer + j*(imax+1) + i + 1;
2464 parent_element = number_of_bulk_elements + (j-1)*imax + i;
2465 fprintf(ptFile,"%d %d %d 0 202 %d %d\n",
2466 n+m, idnr, parent_element + 1, nodes_of_element[0] + 1, nodes_of_element[1] + 1);
2467 }
2468 /* east */
2469 ++idnr;
2470 i=imax;
2471 printf("| western (x=max) boundary for sidebulk; idnr=%d\n", idnr);
2472 for (j=0;j<jmax;++j){ /* loop over all j-values */
2473 ++m;
2474 nodes_of_element[0] = (kmax-1)*nodes_in_one_layer + j*(imax+1) + i;
2475 nodes_of_element[1] = (kmax-1)*nodes_in_one_layer + j*(imax+1) + i+1;
2476 parent_element = number_of_bulk_elements + j*imax + i-1;
2477 fprintf(ptFile,"%d %d %d 0 202 %d %d\n",
2478 n+m, idnr, parent_element + 1, nodes_of_element[0] + 1, nodes_of_element[1] + 1);
2479 }
2480 }
2481 fclose(ptFile);
2482 if (n+m != number_of_boundary_elements + sidebulk*2*(imax+jmax)){
2483 printf("\a %d=%d + %d written boundary elements does not match %d= %d + %d calculated elements in grid\n",
2484 n+m, n,m, number_of_boundary_elements+sidebulk*2*(imax+jmax), number_of_boundary_elements, sidebulk*2*(imax+jmax));
2485 return;
2486 }else{
2487 printf("| %d number of boundary elements for bulk written.\n", n);
2488 printf("| %d number of boundary elements for sidebulk written.\n", m);
2489 printf("| succeeded in writting boundary element file %s.\n",filename);
2490 printf("---------------------------------------------------------------\n");
2491 }
2492
2493 free(iced);free(delta_z);free(bottom);
2494 return;
2495 }
2496
2497
2498 /****************************************************************
2499 Get parameters from log file of SICOPOLIS run
2500 *****************************************************************/
FC_FUNC_(readlog_c,READLOG_C)2501 void STDCALLBULL FC_FUNC_(readlog_c,READLOG_C) (FC_CHAR_PTR(runname,l1), /*name of run*/
2502 int *imax, /* grid steps in xi-direction */
2503 int *jmax, /* grid steps in eta-direction */
2504 int *kcmax, /* grid steps in z-direction cold ice layer */
2505 int *ktmax,/* grid steps in z-direction temperate ice layer */
2506 int *krmax,/* grid steps in z-direction in bedrock */
2507 float *deform, /* parameter for vertical scaling */
2508 float *deltaX, /* horizontal grod spacing */
2509 int *gotit){
2510 /* variable declaration */
2511 int gotparameter=0, i=0, inputint;
2512 char filename[80], inputstring[100], parameter[6], chardummy[1];
2513 float inputfloat;
2514 FILE *ptFile;
2515
2516
2517 /* compose filename */
2518 sprintf(filename,"%s%s", runname, ".log");
2519
2520 /* print out little summary */
2521 printf("---------------------------------------------------------------\n");
2522 printf("| Reading in parameters from log file\n");
2523 printf("---------------------------------------------------------------\n");
2524
2525 if((ptFile = fopen(filename, "r")) == NULL){
2526 printf("Cannot open file %s\n", filename);
2527 *gotit = 0;
2528 return;
2529 }
2530
2531 while((fgets(inputstring, LINESIZE - 1, ptFile) != NULL) && (gotparameter <7)){
2532 ++i;
2533 /* printf("%s\n", inputstring); */
2534 if (sscanf(inputstring, "%s %c %i3", parameter, chardummy, &inputint) == 3){
2535 /* printf("%d: %s, %c, %f\n", i, parameter, chardummy, inputfloat); */
2536 if (strcmp(parameter, "imax") == 0){
2537 *imax = inputint;
2538 gotparameter++;
2539 }
2540 else if (strcmp(parameter, "jmax") == 0){
2541 *jmax = inputint;
2542 gotparameter++;
2543 }
2544 else if (strcmp(parameter, "kcmax") == 0){
2545 *kcmax = inputint;
2546 gotparameter++;
2547 }
2548 else if (strcmp(parameter, "ktmax") == 0){
2549 *ktmax = inputint;
2550 gotparameter++;
2551 }
2552 else if (strcmp(parameter, "krmax") == 0){
2553 *ktmax = inputint;
2554 gotparameter++;
2555 }
2556 }
2557 if (sscanf(inputstring, "%s %c %f", parameter, chardummy, &inputfloat) == 3){
2558 if (strcmp(parameter, "a") == 0){
2559 *deform = inputfloat;
2560 gotparameter++;
2561 }
2562 else if (strcmp(parameter, "dx") == 0){
2563 *deltaX = inputfloat;
2564 gotparameter++;
2565 }
2566 }
2567 }
2568 printf("| imax/jmax/kcmax/ktmax=%d/%d/%d/%d\n",*imax, *jmax, *kcmax, *ktmax);
2569 printf("| deform = %f, dx = %f\n", *deform, *deltaX);
2570 printf("---------------------------------------------------------------\n");
2571 fclose(ptFile);
2572 *gotit=1;
2573 return;
2574 }
2575