TheBoussinesqModel  3.2.1
 All Data Structures Files Functions Variables Typedefs Macros Pages
rw_maps.c
Go to the documentation of this file.
1 
2 /* STATEMENT:
3 
4 ASCII-GEOtop LIBRARIES
5 
6 Copyright, 2008 Stefano Endrizzi, Riccardo Rigon
7 
8  LICENSE:
9 
10  This file is part of ASCII-GEOtop LIBRARIES
11  ASCII-GEOtop LIBRARIES is a free software: you can redistribute it and/or modify
12  it under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  This program is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  GNU General Public License for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with this program. If not, see <http://www.gnu.org/licenses/>.*/
23 
24 
25 #include "turtle.h"
26 #include "tensor3D.h" /* line added by Emanuele Cordano on 1 September 2009 */
27 #include "t_utilities.h" /* line added by Emanuele Cordano on 1 September 2009 */
28 #include "rw_maps.h"
29 #include "import_ascii.h"
30 #include "write_ascii.h"
31 #include "extensions.h"
32 
33 //----------------------------------------------------------------------------------------------
34 //----------------------------------------------------------------------------------------------
35 
36 SHORTMATRIX *copyshort_doublematrix(DOUBLEMATRIX *M){
37 
38  SHORTMATRIX *S;
39  long r, c;
40 
41  S=new_shortmatrix(M->nrh,M->nch);
42  for(r=1;r<=M->nrh;r++){
43  for(c=1;c<=M->nch;c++){
44  S->co[r][c]=(short)M->co[r][c];
45  }
46  }
47 
48  return(S);
49 
50 }
51 
52 //----------------------------------------------------------------------------------------------
53 //----------------------------------------------------------------------------------------------
54 
55 LONGMATRIX *copylong_doublematrix(DOUBLEMATRIX *M){
56 
57  LONGMATRIX *L;
58  long r, c;
59 
60  L=new_longmatrix(M->nrh,M->nch);
61  for(r=1;r<=M->nrh;r++){
62  for(c=1;c<=M->nch;c++){
63  L->co[r][c]=(long)M->co[r][c];
64  }
65  }
66 
67  return(L);
68 
69 }
70 
71 //----------------------------------------------------------------------------------------------
72 //----------------------------------------------------------------------------------------------
73 
74 DOUBLEMATRIX *copydouble_shortmatrix(SHORTMATRIX *S){
75 
76  DOUBLEMATRIX *M;
77  long r, c;
78 
79  M=new_doublematrix(S->nrh,S->nch);
80  for(r=1;r<=S->nrh;r++){
81  for(c=1;c<=S->nch;c++){
82  M->co[r][c]=(double)S->co[r][c];
83  }
84  }
85 
86  return(M);
87 
88 }
89 
90 
91 //----------------------------------------------------------------------------------------------
92 //----------------------------------------------------------------------------------------------
93 
94 DOUBLEMATRIX *copydouble_longmatrix(LONGMATRIX *L){
95 
96  DOUBLEMATRIX *M;
97  long r, c;
98 
99  M=new_doublematrix(L->nrh,L->nch);
100  for(r=1;r<=L->nrh;r++){
101  for(c=1;c<=L->nch;c++){
102  M->co[r][c]=(double)L->co[r][c];
103  }
104  }
105 
106  return(M);
107 
108 }
109 
110 //----------------------------------------------------------------------------------------------
111 //----------------------------------------------------------------------------------------------
112 
113 DOUBLEMATRIX *copydoublematrix_const(double c0, DOUBLEMATRIX *Mref, double NOVALUE){
114 
115  DOUBLEMATRIX *M;
116  long r, c;
117 
118  M=new_doublematrix(Mref->nrh,Mref->nch);
119  for(r=1;r<=M->nrh;r++){
120  for(c=1;c<=M->nch;c++){
121  if(Mref->co[r][c]==NOVALUE){
122  M->co[r][c]=NOVALUE;
123  }else{
124  M->co[r][c]=c0;
125  }
126  }
127  }
128 
129  return(M);
130 
131 }
132 
133 //----------------------------------------------------------------------------------------------
134 //----------------------------------------------------------------------------------------------
135 
136 DOUBLEMATRIX *multiplydoublematrix(double f, DOUBLEMATRIX *Mref, double NOVALUE){
137 
138  DOUBLEMATRIX *M;
139  long r, c;
140 
141  M=new_doublematrix(Mref->nrh,Mref->nch);
142  for(r=1;r<=M->nrh;r++){
143  for(c=1;c<=M->nch;c++){
144  if(Mref->co[r][c]==NOVALUE){
145  M->co[r][c]=NOVALUE;
146  }else{
147  M->co[r][c]=f*Mref->co[r][c];
148  }
149  }
150  }
151 
152  return(M);
153 
154 }
155 
156 //----------------------------------------------------------------------------------------------
157 //----------------------------------------------------------------------------------------------
158 
159 void build_doubletensor(DOUBLETENSOR *T, DOUBLEMATRIX *M, long l){
160 
161  long r,c;
162 
163  if(l<=0 || l>T->ndh) t_error("Invalid doubletensor construction");
164  if(T->nrh!=M->nrh) t_error("Invalid doubletensor construction");
165  if(T->nch!=M->nch) t_error("Invalid doubletensor construction");
166 
167  for(r=1;r<=M->nrh;r++){
168  for(c=1;c<=M->nch;c++){
169  T->co[l][r][c]=M->co[r][c];
170  }
171  }
172 
173 }
174 
175 //----------------------------------------------------------------------------------------------
176 //----------------------------------------------------------------------------------------------
177 
178 DOUBLEMATRIX *extract_doublematrix(DOUBLETENSOR *T, long l){
179 
180  long r,c;
181  DOUBLEMATRIX *M;
182 
183  if(l<=0 || l>T->ndh) t_error("Invalid doubletensor extraction");
184 
185  M=new_doublematrix(T->nrh,T->nch);
186 
187  for(r=1;r<=M->nrh;r++){
188  for(c=1;c<=M->nch;c++){
189  M->co[r][c]=T->co[l][r][c];
190  }
191  }
192 
193  return(M);
194 
195 }
196 
197 
198 
199 //----------------------------------------------------------------------------------------------
200 //----------------------------------------------------------------------------------------------
201 
202 DOUBLEMATRIX *extract_fromtensor(DOUBLETENSOR *T, long l){
203 
204  long r, c;
205  DOUBLEMATRIX *M;
206 
207  if(l<1 || l>T->ndh) t_error("cannot extract a matrix from a tensor");
208 
209  M=new_doublematrix(T->nrh,T->nch);
210  for(r=1;r<=T->nrh;r++){
211  for(c=1;c<=T->nch;c++){
212  M->co[r][c]=T->co[l][r][c];
213  }
214  }
215 
216  return(M);
217 
218 }
219 
220 //----------------------------------------------------------------------------------------------
221 //----------------------------------------------------------------------------------------------
222 
223 DOUBLETENSOR *build_frommatrix(DOUBLEMATRIX *M, long l, long lmax){
224 
225  long r,c;
226  DOUBLETENSOR *T;
227 
228  if(lmax<1) t_error("cannot build a tensor from a matrix");
229  if(l>lmax) t_error("cannot build a tensor from a matrix");
230 
231  T=new_doubletensor(lmax, M->nrh, M->nch);
232  initialize_doubletensor(T,0.0);
233  for(r=1;r<=M->nrh;r++){
234  for(c=1;c<=M->nch;c++){
235  T->co[l][r][c]=M->co[r][c];
236  }
237  }
238 
239  return(T);
240 
241 }
242 
243 
244 //----------------------------------------------------------------------------------------------
245 //----------------------------------------------------------------------------------------------
246 
247 void write_frommatrix(long l, DOUBLEMATRIX *M, DOUBLETENSOR *T){
248 
249  long r,c;
250 
251  if(l<0 || l>T->ndh) t_error("cannot write a tensor from a matrix");
252  if(M->nrh!=T->nrh) t_error("cannot write a tensor from a matrix");
253  if(M->nch!=T->nch) t_error("cannot write a tensor from a matrix");
254 
255  for(r=1;r<=M->nrh;r++){
256  for(c=1;c<=M->nch;c++){
257  T->co[l][r][c]=M->co[r][c];
258  }
259  }
260 
261 }
262 
263 
264 //----------------------------------------------------------------------------------------------
265 //----------------------------------------------------------------------------------------------
266 
267 void fmultiplydoublematrix(DOUBLEMATRIX *destination, DOUBLEMATRIX *origin, double f, double novalue){
268 
269  long r,c;
270 
271  if(origin->nrh!=destination->nrh) t_error("cannot copy matrix");
272  if(origin->nch!=destination->nch) t_error("cannot copy matrix");
273  for(r=1;r<=origin->nrh;r++){
274  for(c=1;c<=origin->nch;c++){
275  if(origin->co[r][c]!=novalue){
276  destination->co[r][c]=f*origin->co[r][c];
277  }else{
278  destination->co[r][c]=novalue;
279  }
280  }
281  }
282 }
283 
284 //----------------------------------------------------------------------------------------------
285 //----------------------------------------------------------------------------------------------
286 
287 void assignnovalue(DOUBLEMATRIX *destination, DOUBLEMATRIX *origin, double novalue){
288 
289  long r,c;
290 
291  if(origin->nrh!=destination->nrh) t_error("cannot assign novalue");
292  if(origin->nch!=destination->nch) t_error("cannot assign novalue");
293  for(r=1;r<=origin->nrh;r++){
294  for(c=1;c<=origin->nch;c++){
295  if(origin->co[r][c]==novalue) destination->co[r][c]=novalue;
296  }
297  }
298 }
299 
300 
301 //----------------------------------------------------------------------------------------------
302 //----------------------------------------------------------------------------------------------
303 
304 void write_suffix(char *suffix, long i, short start){
305 
306  short m, c, d, u;
307 
308  if(i>=0 && i<=9){
309  m=0;
310  c=0;
311  d=0;
312  u=i;
313  }else if(i<=99){
314  m=0;
315  c=0;
316  d=(short)(i/10.0);
317  u=i-10.0*d;
318  }else if(i<=999){
319  m=0;
320  c=(short)(i/100.0);
321  d=(short)((i-100.0*c)/10.0);
322  u=i-100.0*c-10*d;
323  }else if(i<=9999){
324  m=(short)(i/1000.0);
325  c=(short)((i-1000.0*m)/100.0);
326  d=(short)((i-1000.0*m-100.0*c)/10.0);
327  u=i-1000*m-100.0*c-10*d;
328  }else{
329  t_error("Number too high");
330  }
331 
332  m+=48;
333  c+=48;
334  d+=48;
335  u+=48;
336 
337  suffix[start]=m;
338  suffix[start+1]=c;
339  suffix[start+2]=d;
340  suffix[start+3]=u;
341 
342 }
343 
344 //----------------------------------------------------------------------------------------------
345 //----------------------------------------------------------------------------------------------
346 
347 char *namefile_i(char *name, long i){
348 
349  char SSSS[ ]={"SSSS"};
350  char *name_out;
351 
352  write_suffix(SSSS, i, 0);
353 
354  name_out=join_strings(name,SSSS);
355  name_out=join_strings(name_out,textfile);
356 
357  return(name_out);
358 
359 }
360 
361 char *namefile_i_we(char *name, long i){
362 
363  char SSSS[ ]={"LSSSS"};
364  char *name_out;
365 
366  write_suffix(SSSS, i, 1);
367 
368  name_out=join_strings(name,SSSS);
369 
370  return(name_out);
371 
372 }
373 
374 char *namefile_i_we2(char *name, long i){
375 
376  char SSSS[ ]={"SSSS"};
377  char *name_out;
378 
379  write_suffix(SSSS, i, 0);
380 
381  name_out=join_strings(name,SSSS);
382 
383  return(name_out);
384 
385 }
386 
387 //----------------------------------------------------------------------------------------------
388 //----------------------------------------------------------------------------------------------
389 
390 short existing_file(char *name){
391 
392  //if the file exists gives 1 (fluidturtle), 2(grassascii), 3(esriascii), 0 if the file doesn't exist
393 
394  short a=0;
395  FILE *f;
396  char *ft, *grass, *esri;
397 
398  ft=join_strings(name,ascii_ft);
399  grass=join_strings(name,ascii_grass);
400  esri=join_strings(name,ascii_esri);
401  if( (f=fopen(ft,"r"))!=NULL ){
402  a=1;
403  fclose(f);
404  }else if( (f=fopen(grass,"r"))!=NULL ){
405  a=2;
406  fclose(f);
407  }else if( (f=fopen(esri,"r"))!=NULL ){
408  a=3;
409  fclose(f);
410  } else {
411  a=-1; /* no file in an acceptable format exists (modified by Emanuele Cordano on August 8, 2009*/
412  }
413  free(ft);free(grass);free(esri);
414  return(a);
415 
416 }
417 
418 short existing_file_text(char *name){
419 
420  //if the file exists gives 1, 0 if the file doesn't exist
421 
422  short a=0;
423  FILE *f;
424 
425  if( (f=fopen(join_strings(name,textfile),"r"))!=NULL ){
426  a=1;
427  fclose(f);
428  }
429  return(a);
430 
431 }
432 //----------------------------------------------------------------------------------------------
433 //----------------------------------------------------------------------------------------------
434 
435 DOUBLEMATRIX *read_map(short a, char *filename, DOUBLEMATRIX *Mref, T_INIT *UVref){
436 
437 // a=0 non usa Mref, UVref output
438 // a=1 non esegue controllo non values, Mref e UVref input
439 // a=2 esegue controllo novalues, Mref e UVref input
440 
441  DOUBLEMATRIX *M;
442  FILE *f;
443  long index, r, c, nr, nc;
444  double *header, *m;
445  double Dxmap, Dymap, X0map, Y0map, novalue;
446  char *filename2;
447  char *ft,*grass,*esri;
448 
449  ft=join_strings(filename,ascii_ft);
450  grass=join_strings(filename,ascii_grass);
451  esri=join_strings(filename,ascii_esri);
452 
453  T_INIT *UV;
454 
455  if (existing_file(filename)==-1) return NULL;
456 
457  if (existing_file(filename)==1){
458  filename2=ft;
459  if(a>0){
460  UV=(T_INIT *)malloc(sizeof(T_INIT));
461  if(!UV) t_error("UV was not allocated in read_map");
462  }
463  f=t_fopen(ft,"r");
464  index=read_index(f,PRINT);
465  if(a==0){
466  UVref->U=read_doublearray(f,PRINT);
467  UVref->V=read_doublearray(f,PRINT);
468  }else{
469  UV->U=read_doublearray(f,PRINT);
470  UV->V=read_doublearray(f,PRINT);
471  }
472  M=read_doublematrix(f,"a",PRINT);
473  t_fclose(f);
474 
475  if(a>0){
476  //Check header
477  if(UV->U->co[2]!=UVref->U->co[2]){
478  printf("Dx in %s file is different from Dx in DTM file! \n",filename2);
479  t_error("Inconsistent map");
480  }
481  if(UV->U->co[1]!=UVref->U->co[1]){
482  printf("Dy in %s file is different from Dy in DTM file! \n",filename2);
483  t_error("Inconsistent map");
484  }
485  if(UV->U->co[4]!=UVref->U->co[4]){
486  printf("X0 in %s file is different from X0 in DTM file! \n",filename2);
487  t_error("Inconsistent map");
488  }
489  if(UV->U->co[3]!=UVref->U->co[3]){
490  printf("Y0 in %s file is different from Y0 in DTM file! \n",filename2);
491  t_error("Inconsistent map");
492  }
493  if(M->nrh!=Mref->nrh){
494  printf("Number of rows in %s file is not consistent with DTM file! \n",filename2);
495  t_error("Inconsistent map");
496  }
497  if(M->nch!=Mref->nch){
498  printf("Number of columns in %s file is not consistent with DTM file! \n",filename2);
499  t_error("Inconsistent map");
500  }
501 
502  //The novalue is imposed egual to that of DTM file:*/
503  for(r=1;r<=M->nrh;r++){
504  for(c=1;c<=M->nch;c++){
505  if(M->co[r][c]==UV->V->co[2]) M->co[r][c]=UVref->V->co[2];
506  if(a>1){
507  if(M->co[r][c]==UVref->V->co[2] && Mref->co[r][c]!=UVref->V->co[2]){
508  printf("Novalues not consistent in %s file",filename2);
509  printf("\nr:%ld c:%ld Mref:%f M:%f\n",r,c,Mref->co[r][c], M->co[r][c]);
510  t_error("Inconsistent map");
511  }
512  if(M->co[r][c]!=UVref->V->co[2] && Mref->co[r][c]==UVref->V->co[2]) M->co[r][c]=UVref->V->co[2];
513  }else{
514  if(M->co[r][c]==UVref->V->co[2] && Mref->co[r][c]!=UVref->V->co[2]){
515  printf("Novalues not consistent in %s file",filename2);
516  printf("\nr:%ld c:%ld Mref:%f M:%f\n",r,c,Mref->co[r][c], M->co[r][c]);
517  stop_execution();
518  }
519  }
520  }
521  }
522 
523  free_doublevector(UV->U);
524  free_doublevector(UV->V);
525  free(UV);
526  }
527 
528  }else{
529  if(a>0){
530  novalue=UVref->V->co[2];
531  }else{
532  novalue=-9999.0;
533  }
534  header = (double *) malloc(6*sizeof(double));
535 
536  if(existing_file(filename)==2){
537  m=read_grassascii(header, novalue, filename, 30);
538  nr=(long)header[4];
539  nc=(long)header[5];
540  Dxmap=(header[2]-header[3])/((long)header[5]);
541  Dymap=(header[0]-header[1])/((long)header[4]);
542  X0map=header[3];
543  Y0map=header[1];
544  filename2=grass;
545  }else if(existing_file(filename)==3){
546  m=read_esriascii(header, novalue, filename, 30);
547  nr=(long)header[1];
548  nc=(long)header[0];
549  Dxmap=header[4];
550  Dymap=header[4];
551  X0map=header[2];
552  Y0map=header[3];
553  filename2=esri;
554  }else{
555  printf("The file %s doesn't exist\n",filename);
556  t_error("Fatal error");
557  }
558 
559  free(header);
560  if(a>0){
561  //Check header
562  if(Dxmap!=UVref->U->co[2]){
563  printf("Dx in %s file is different from Dx in DTM file! \n",filename2);
564  t_error("Inconsistent map");
565  }
566  if(Dymap!=UVref->U->co[1]){
567  printf("Dy:%f in %s file is different from Dy:%f in DTM file! \n",Dymap,filename2,UVref->U->co[1]);
568  t_error("Inconsistent map");
569  }
570  if(X0map!=UVref->U->co[4]){
571  printf("X0 in %s file is different from X0 in DTM file! \n",filename2);
572  t_error("Inconsistent map");
573  }
574  if(Y0map!=UVref->U->co[3]){
575  printf("Y0 in %s file is different from Y0 in DTM file! \n",filename2);
576  t_error("Inconsistent map");
577  }
578  if(nr!=Mref->nrh){
579  printf("Number of rows in %s file (%ld) is not consistent with DTM file (%ld)! \n",filename2,nr,Mref->nrh);
580  t_error("Inconsistent map");
581  }
582  if(nc!=Mref->nch){
583  printf("Number of columns in %s file is not consistent with DTM file! \n",filename2);
584  t_error("Inconsistent map");
585  }
586 
587  }else{
588  UVref->U=new_doublevector(4);
589  UVref->V=new_doublevector(2);
590  UVref->U->co[2]=Dxmap;
591  UVref->U->co[1]=Dymap;
592  UVref->U->co[4]=X0map;
593  UVref->U->co[3]=Y0map;
594  UVref->V->co[1]=-1;
595  UVref->V->co[2]=novalue;
596  }
597 
598  //assign values and check novalues
599  M=new_doublematrix(nr,nc);
600  for(r=1;r<=nr;r++){
601  for(c=1;c<=nc;c++){
602  M->co[r][c]=m[(r-1)*nc+c-1];
603  if(a>1){
604  if (M->co[r][c]==UVref->V->co[2] && Mref->co[r][c]!=UVref->V->co[2]){
605  printf("Novalues not consistent in %s file",filename2);
606  printf("\nr:%ld c:%ld Mref:%f M:%f\n",r,c,Mref->co[r][c], M->co[r][c]);
607  t_error("Inconsistent map");
608  }
609  if(M->co[r][c]!=UVref->V->co[2] && Mref->co[r][c]==UVref->V->co[2]) M->co[r][c]=UVref->V->co[2];
610  }
611  }
612  }
613 
614  free(m);
615  }
616  free(ft);free(grass);free(esri);
617  return(M);
618 
619 }
620 
621 
622 //----------------------------------------------------------------------------------------------
623 //----------------------------------------------------------------------------------------------
624 
625 DOUBLEMATRIX *read_mapseries(long i, char *filename, DOUBLEMATRIX *Mref, T_INIT *UVref){
626 
627  char SSSS[ ]={"LSSSS"};
628  char *name;
629  DOUBLEMATRIX *M;
630 
631  write_suffix(SSSS, i, 1);
632  name=join_strings(filename,SSSS);
633  M=read_map(2, name, Mref, UVref);
634 
635  return(M);
636 
637 }
638 
639 //----------------------------------------------------------------------------------------------
640 //----------------------------------------------------------------------------------------------
641 
642 DOUBLETENSOR *read_tensor(long nl, char *filename, DOUBLEMATRIX *Mref, T_INIT *UVref){
643  long l;
644  DOUBLEMATRIX *M;
645  DOUBLETENSOR *T;
646 
647  if(nl<1) t_error("The dimension of the tensor must be greater than or equal to 1");
648 
649  for(l=1;l<=nl;l++){
650  M=read_mapseries(l, filename, Mref, UVref);
651  if(l==1){
652  T=build_frommatrix(M, l, nl);
653  }else{
654  write_frommatrix(l, M, T);
655  }
656  free_doublematrix(M);
657  }
658 
659  return(T);
660 }
661 
662 
663 
664 //----------------------------------------------------------------------------------------------
665 //----------------------------------------------------------------------------------------------
666 
667 DOUBLETENSOR *read_maptensor(long i, long lmax, char *filename, DOUBLEMATRIX *Mref, T_INIT *UVref){
668 
669  char SSSSLLLLL[ ]={"SSSSLLLLL"};
670  char *name;
671  DOUBLETENSOR *T;
672  DOUBLEMATRIX *M;
673  long l;
674 
675  if(i<0 || lmax<1) t_error("cannot read a tensor with null or negative columns");
676 
677  write_suffix(SSSSLLLLL, i, 0);
678 
679  for(l=1;l<=lmax;l++){
680  write_suffix(SSSSLLLLL, l, 5);
681  name=join_strings(filename,SSSSLLLLL);
682  M=read_map(2, name, Mref, UVref);
683  if(l==1){
684  T=build_frommatrix(M, l, lmax);
685  }else{
686  write_frommatrix(l, M, T);
687  }
688  }
689 
690  return(T);
691 
692 }
693 
694 //----------------------------------------------------------------------------------------------
695 //----------------------------------------------------------------------------------------------
696 
697 void write_map(char *filename, short type, short format, DOUBLEMATRIX *M, T_INIT *UV){
698 
699 // type=0 floating point
700 // type=1 integer
701 
702 // format=1 fluidturtle
703 // format=2 grassascii
704 // format=3 esriascii
705 
706  if(format==1){
707  write_fluidturtle(filename, type, M, UV);
708  }else if(format==2){
709  write_grassascii(filename, type, M, UV);
710  }else if(format==3){
711  write_esriascii(filename, type, M, UV);
712  }
713 
714 }
715 
716 //----------------------------------------------------------------------------------------------
717 //----------------------------------------------------------------------------------------------
718 
719 void write_mapseries(long i, char *filename, short type, short format, DOUBLEMATRIX *M, T_INIT *UV){
720 
721  char SSSS[ ]={"SSSS"};
722 
723  write_suffix(SSSS, i, 0);
724  write_map(join_strings(filename, SSSS), type, format, M, UV);
725 
726 }
727 
728 //----------------------------------------------------------------------------------------------
729 //----------------------------------------------------------------------------------------------
730 
731 void write_tensorseries(short a, long l, long i, char *filename, short type, short format, DOUBLETENSOR *T, T_INIT *UV){
732 
733 // a=0 non include "l" nel suffisso
734 // a=1 include "l" nel suffisso
735 // l:layer
736 // i:temporal step
737 
738  char SSSSLLLLL[ ]={"SSSSLLLLL"};
739  char SSSS[ ]={"SSSS"};
740  char *name;
741  long r, c;
742  DOUBLEMATRIX *M;
743 
744  if(a==0){
745  write_suffix(SSSS, i, 0);
746  name=join_strings(filename,SSSS);
747  }else if(a==1){
748  write_suffix(SSSSLLLLL, i, 0);
749  write_suffix(SSSSLLLLL, l, 5);
750  name=join_strings(filename,SSSSLLLLL);
751  }
752 
753  M=new_doublematrix(T->nrh,T->nch);
754  for(r=1;r<=T->nrh;r++){
755  for(c=1;c<=T->nch;c++){
756  M->co[r][c]=T->co[l][r][c];
757  }
758  }
759 
760  write_map(name, type, format, M, UV);
761 
762  free_doublematrix(M);
763 
764 }
765 
766 void write_tensorseries_bis(short a, long l, long i, char *filename, short type, short format, DOUBLETENSOR *T, T_INIT *UV){
767 
768 // a=0 non include "l" nel suffisso
769 // a=1 include "l" nel suffisso
770 // l:layer
771 // i:temporal step
772 
773  char SSSSLLLLL[ ]={"LLLLLNNNNN"};
774  char SSSS[ ]={"NNNN"};
775  char *name;
776  long r, c;
777  DOUBLEMATRIX *M;
778 
779  if(a==0){
780  write_suffix(SSSS, i, 0);
781  name=join_strings(filename,SSSS);
782  }else if(a==1){
783  write_suffix(SSSSLLLLL, l, 1);
784  write_suffix(SSSSLLLLL, i, 6);
785  name=join_strings(filename,SSSSLLLLL);
786  }
787 
788  M=new_doublematrix(T->nrh,T->nch);
789  for(r=1;r<=T->nrh;r++){
790  for(c=1;c<=T->nch;c++){
791  M->co[r][c]=T->co[l][r][c];
792  }
793  }
794 
795  write_map(name, type, format, M, UV);
796 
797  free_doublematrix(M);
798 
799 }
800 
801 //----------------------------------------------------------------------------------------------
802 //----------------------------------------------------------------------------------------------
803 
804 void write_tensorseries2(long i, char *filename, short type, short format, DOUBLETENSOR *T, T_INIT *UV){
805 
806  long l;
807  for(l=1;l<=T->ndh;l++){
808  write_tensorseries_bis(1, l, i, filename, type, format, T, UV);
809  }
810 }
811 
812 //----------------------------------------------------------------------------------------------
813 //----------------------------------------------------------------------------------------------
814