GRIDMAN
grid managment library
grid1.f
Go to the documentation of this file.
1 C> @file grid1.f
2 C> Basic methods of the data-type GRIDMAN_GRID, see gridman.f
3 C GRIDMAN, grid managment library. Author: Vladislav Kotov, v.kotov@fz-juelich.de
4 
5 ! Copyright (c) 2017 Forschungszentrum Juelich GmbH
6 ! Vladislav Kotov
7 !
8 ! This file is part of GRIDMAN.
9 !
10 ! GRIDMAN is free software: you can redistribute it and/or modify
11 ! it under the terms of the GNU General Public License as published by
12 ! the Free Software Foundation, either version 3 of the License, or
13 ! (at your option) any later version.
14 !
15 ! GRIDMAN is distributed in the hope that it will be useful,
16 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 ! GNU General Public License for more details.
19 !
20 ! You should have received a copy of the GNU General Public License
21 ! along with GRIDMAN. If not, see <http://www.gnu.org/licenses/>.
22 
23 C *****************************************************************************************
24 C> Allocate GRIDMAN_GRID object
25 C *****************************************************************************************
26 C>
27 C> WARNING: if GRID exists it will be overwritten!
28  SUBROUTINE gridman_grid_allocate(GRID,TYPE,NEDGES,NPOINTS,NCELLS,
29  s ierr,nedgeindex,ncellindex)
32  IMPLICIT NONE
33  INTRINSIC PRESENT
34 C> Resulting grid object
35  TYPE(gridman_grid) :: GRID
36 C> Grid type
37  INTEGER,INTENT(IN) :: TYPE
38 C> Number of cells
39  INTEGER(GRIDMAN_SP),INTENT(IN) :: NCELLS
40 C> Number of edges
41  INTEGER(GRIDMAN_SP),INTENT(IN) :: NEDGES
42 C> Number of points
43  INTEGER(GRIDMAN_SP),INTENT(IN) :: NPOINTS
44 C> Error code
45  INTEGER,INTENT(OUT) :: IERR
46 C> Number of groups of indices for edges
47 C>
48 C> If NEDGEINDEX<=0 then egde indices are not allocated
49  INTEGER,INTENT(IN),OPTIONAL :: NEDGEINDEX
50 C> Number of groups of indices for cells
51 C>
52 C> If NEDGEINDEX<=0 then cell indices are not allocated
53  INTEGER,INTENT(IN),OPTIONAL :: NCELLINDEX
54 
55  INTEGER :: ST,IERR0
56 
57  IF(gridman_dbg)
58  w WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_ALLOCATE"
59 
60  ierr=0
61 
62  CALL gridman_grid_deallocate(grid,ierr)
63  IF(ierr.NE.0) GOTO 1000
64 
65  grid%TYPE=TYPE
66  CALL dimensions(TYPE,GRID%PDIM,GRID%EDIM,IERR)
67  IF(ierr.GT.0) GOTO 1000
68 
69  IF(nedges.LT.1) THEN
70  WRITE(gridman_unit,*)
71  w "ERROR in GRIDMAN_GRID_ALLOCATE: incorrect number of edges"
72  WRITE(gridman_unit,*) " NEDGES ",nedges
73  ierr=100
74  RETURN
75  END IF
76  grid%NEDGES=nedges
77  ALLOCATE(grid%CELLS(2,nedges),
78  a grid%POINTS(grid%EDIM,nedges),stat=st)
79  IF(st.NE.0) THEN
80  WRITE(gridman_unit,*)
81  w "ERROR in GRIDMAN_GRID_ALLOCATE: cannot allocate edge arrays"
82  WRITE(gridman_unit,*) " NEDGES ",nedges
83  ierr=200
84  RETURN
85  END IF
86 
87  IF(npoints.LT.1) THEN
88  WRITE(gridman_unit,*)
89  w "ERROR in GRIDMAN_GRID_ALLOCATE: incorrect number of points"
90  WRITE(gridman_unit,*) " NPOINTS ",npoints
91  ierr=100
92  RETURN
93  END IF
94  grid%NPOINTS=npoints
95  ALLOCATE(grid%X(grid%PDIM,npoints),stat=st)
96  IF(st.NE.0) THEN
97  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ALLOCATE: ",
98  w "cannot allocate coordinate arrays"
99  WRITE(gridman_unit,*) " NPOINTS ",npoints
100  ierr=200
101  RETURN
102  END IF
103  grid%NCELLS=ncells
104 
105 C INDICES
106  IF(PRESENT(nedgeindex)) THEN
107  IF(nedgeindex.LT.1) THEN
108  grid%NEDGEINDEX=0
109  ELSE
110  grid%NEDGEINDEX=nedgeindex
111  ALLOCATE(grid%EDGEINDEX(nedgeindex),stat=st)
112  IF(st.NE.0) THEN
113  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ALLOCATE: ",
114  w "cannot allocate edge indices"
115  WRITE(gridman_unit,*) " NEDGEINDEX ",nedgeindex
116  ierr=200
117  END IF
118  END IF
119  ELSE
120  grid%NEDGEINDEX=0
121  END IF
122 
123  IF(PRESENT(ncellindex)) THEN
124  IF(ncellindex.LT.1) THEN
125  grid%NCELLINDEX=0
126  ELSE
127  grid%NCELLINDEX=ncellindex
128  ALLOCATE(grid%CELLINDEX(ncellindex),stat=st)
129  IF(st.NE.0) THEN
130  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_ALLOCATE: ",
131  w "cannot allocate edge indices"
132  WRITE(gridman_unit,*) " NCELLINDEX ",ncellindex
133  ierr=200
134  END IF
135  END IF
136  ELSE
137  grid%NCELLINDEX=0
138  END IF
139 
140 C DEFAULTS
141  grid%UNIT2SI=1.0
142  grid%UNITS='METER'
143  grid%DESCRIPTION='Created by GRIDMAN_GRID_ALLOCATE'
144 
145  IF(gridman_dbg)
146  w WRITE(gridman_unit,*) "GRIDMAN_GRID_ALLOCATE finished"
147 
148  RETURN
149 
150 1000 WRITE(gridman_unit,*) "GRIDMAN_GRID_ALLOCATE is terminated"
151  CALL gridman_grid_deallocate(grid,ierr0)
152 
153  CONTAINS
154 C
155 C DIMENSIONS AS FUNCTION OF THE GRID TYPE
156 C
157  SUBROUTINE dimensions(TYPE,PDIM,EDIM,IERR)
158 
159  INTEGER,INTENT(IN) :: TYPE
160  INTEGER,INTENT(OUT) :: PDIM,EDIM,IERR
161 
162  SELECT CASE (type)
163  CASE(2)
164  edim=2
165  pdim=2
166  CASE(3)
167  edim=3
168  pdim=3
169  CASE DEFAULT
170  ierr=100
171  WRITE(gridman_unit,*)
172  w "ERROR in GRIDMAN_GRID_DIMENSIONS: unknown grid type"
173  WRITE(gridman_unit,*) " TYPE ",TYPE
174  WRITE(gridman_unit,*) "Known grid types: 2,3"
175  END SELECT
176  END SUBROUTINE dimensions
177 
178  END SUBROUTINE gridman_grid_allocate
179 
180 C*****************************************************************************************
181 C> Deallocate grid object
182 C*****************************************************************************************
183  SUBROUTINE gridman_grid_deallocate(GRID,IERR)
186  IMPLICIT NONE
187  INTRINSIC SIZE
188 C> Grid object to be deallocated
189  TYPE(gridman_grid) :: GRID
190 C> Error code
191  INTEGER,INTENT(OUT) :: IERR
192 
193  INTEGER :: ST,I
194 
195  IF(gridman_dbg)
196  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_DEALLOCATE"
197 
198  ierr=0
199 
200  IF(ALLOCATED(grid%CELLS)) THEN
201  DEALLOCATE(grid%CELLS,stat=st)
202  IF(st.NE.0) THEN
203  ierr=-200
204  WRITE(gridman_unit,*)
205  w "WARNING from GRIDMAN_GRID_DEALLOCATE: ",
206  w "can not deallocate CELLS"
207  END IF
208  END IF
209 
210  IF(ALLOCATED(grid%POINTS)) THEN
211  DEALLOCATE(grid%POINTS,stat=st)
212  IF(st.NE.0) THEN
213  ierr=-200
214  WRITE(gridman_unit,*)
215  w "WARNING from GRIDMAN_GRID_DEALLOCATE: ",
216  w "can not deallocate POINTS"
217  END IF
218  END IF
219 
220  IF(ALLOCATED(grid%X)) THEN
221  DEALLOCATE(grid%X,stat=st)
222  IF(st.NE.0) THEN
223  ierr=-200
224  WRITE(gridman_unit,*)
225  w "WARNING from GRIDMAN_GRID_DEALLOCATE: ",
226  w "can not deallocate coordinates X"
227  END IF
228  END IF
229 
230  IF(ALLOCATED(grid%EDGEINDEX)) THEN
231  DO i=1,SIZE(grid%EDGEINDEX)
232  CALL gridman_index_deallocate(grid%EDGEINDEX(i),ierr)
233  IF(ierr.NE.0) THEN
234  WRITE(gridman_unit,*)
235  w "WARNING from GRIDMAN_GRID_DEALLOCATE: ",
236  w "can not deallocate EDGEINDEX(I), I ",i
237  END IF
238  END DO
239  DEALLOCATE(grid%EDGEINDEX,stat=st)
240  IF(st.NE.0) THEN
241  ierr=-200
242  WRITE(gridman_unit,*)
243  w "WARNING from GRIDMAN_GRID_DEALLOCATE: ",
244  w "can not deallocate EDGEINDEX"
245  END IF
246  END IF
247 
248  IF(ALLOCATED(grid%CELLINDEX)) THEN
249  DO i=1,SIZE(grid%CELLINDEX)
250  CALL gridman_index_deallocate(grid%CELLINDEX(i),ierr)
251  IF(ierr.NE.0) THEN
252  WRITE(gridman_unit,*)
253  w "WARNING from GRIDMAN_GRID_DEALLOCATE: ",
254  w "can not deallocate CELLINDEX(I), I ",i
255  END IF
256  END DO
257  DEALLOCATE(grid%CELLINDEX,stat=st)
258  IF(st.NE.0) THEN
259  ierr=-200
260  WRITE(gridman_unit,*)
261  w "WARNING from GRIDMAN_GRID_DEALLOCATE: ",
262  w "can not deallocate CELLINDEX"
263  END IF
264  END IF
265 
266  IF(gridman_dbg)
267  w WRITE(gridman_unit,*) "GRIDMAN_GRID_DEALLOCATE finished"
268 
269  END SUBROUTINE gridman_grid_deallocate
270 
271 
272 C*****************************************************************************************
273 C> Check consistency of the grid data
274 C*****************************************************************************************
275 C>
276 C> - RES=0: the grid is correct
277 C> - RES=1: grid is not allocated
278 C> - RES=2: dimensions are incorrect
279 C> - RES=3: wrong array size
280 C> - RES=4: cell index is out of bounds
281 C> - RES=5: duplicate cell index for an edge
282 C> - RES=6: NCELLS and CELLS do not match
283 C> - RES=7: index of point is out of bounds
284 C> - RES=8: incorrect unit scaling factor
285 C> - RES=9: incorrect edge indices
286 C> - RES=10: incorrect cell indices
287 C
288  SUBROUTINE gridman_grid_check(GRID,RES,IERR)
290  u gridman_dbg
292  IMPLICIT NONE
293  INTRINSIC SIZE,tiny,ALLOCATED
294 C> Grid object to be checked
295  TYPE(gridman_grid) :: GRID
296 C> Code of the result
297  INTEGER,INTENT(OUT) :: RES
298 C> Error code
299  INTEGER,INTENT(OUT) :: IERR
300 
301  INTEGER (GRIDMAN_SP) :: N1,N2,ICELL,I1,I2,IE,IEDGE,NCELLS_MAX
302  INTEGER :: ID
303 
304  IF(gridman_dbg)
305  w WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_CHECK"
306 
307  ierr=0
308 
309  IF(.NOT.ALLOCATED(grid%CELLS).OR.
310  . .NOT.ALLOCATED(grid%POINTS).OR.
311  . .NOT.ALLOCATED(grid%X)) THEN
312  res=1
313  WRITE(gridman_unit,*)
314  w " GRIDMAN_GRID_CHECK: grid is not allocated"
315  GOTO 100
316  END IF
317 
318 C CHECK DIMENSIONS
319  IF(grid%NCELLS.LT.0.OR.grid%NEDGES.LT.1.OR. !NUMBER OF CELLS CAN BE EQUAL ZERO!
320  f grid%NPOINTS.LT.1.OR.grid%EDIM.LT.1.OR.grid%PDIM.LT.1) THEN
321  res=2
322  WRITE(gridman_unit,*)
323  w " GRIDMAN_GRID_CHECK: incorrect dimensions"
324  WRITE(gridman_unit,*)
325  w "NCELL, NEDGES, NPOINTS, PDIM, EDIM ",
326  w grid%NCELLS,grid%NEDGES,grid%NPOINTS,
327  w grid%PDIM,grid%EDIM
328  GOTO 100
329  END IF
330 
331 C CHECK FOR MISMATCH BETWEEN DIMENSIONS AND ARRAY INDICES
332  n1=SIZE(grid%CELLS,1)
333  n2=SIZE(grid%CELLS,2)
334  IF(n1.NE.2.OR.n2.LT.grid%NEDGES) THEN
335  res=3
336  WRITE(gridman_unit,*)
337  w " GRIDMAN_GRID_CHECK: array size mismatch, CELLS"
338  WRITE(gridman_unit,*) " NEDGES, SIZE(CELLS) ",grid%NEDGES,n1,n2
339  GOTO 100
340  END IF
341  n1=SIZE(grid%POINTS,1)
342  n2=SIZE(grid%POINTS,2)
343  IF(n1.NE.grid%EDIM.OR.n2.LT.grid%NEDGES) THEN
344  res=3
345  WRITE(gridman_unit,*)
346  w " GRIDMAN_GRID_CHECK: array size mismatch, POINTS"
347  WRITE(gridman_unit,*) " EDIM, NEDGES, SIZE(POINTS) ",
348  w grid%EDIM,grid%NEDGES,n1,n2
349  GOTO 100
350  END IF
351  n1=SIZE(grid%X,1)
352  n2=SIZE(grid%X,2)
353  IF(n1.NE.grid%PDIM.OR.n2.LT.grid%NPOINTS) THEN
354  res=3
355  WRITE(gridman_unit,*)
356  w " GRIDMAN_GRID_CHECK: array size mismatch, X"
357  WRITE(gridman_unit,*) " PDIM, NPOINTS, SIZE(X) ",
358  w grid%PDIM,grid%NPOINTS,n1,n2
359  GOTO 100
360  END IF
361 
362 C CELL INDICES OUT OF BOUNDS OR DUPLICATE CELL INDICES
363  DO iedge=1,grid%NEDGES
364  i1=grid%CELLS(1,iedge)
365  i2=grid%CELLS(2,iedge)
366  IF(i1.GT.grid%NCELLS.OR.i1.LT.-grid%NEDGES.OR.
367  w i2.GT.grid%NCELLS.OR.i2.LT.-grid%NEDGES) THEN
368  res=4
369  WRITE(gridman_unit,*)
370  w " GRIDMAN_GRID_CHECK: cells out of bounds ",
371  w " IEDGE ,CELLS, NCELLS, NEDGES",
372  w iedge,grid%CELLS(:,iedge),grid%NCELLS,grid%NEDGES
373  GOTO 100
374  END IF
375  IF(i2.EQ.i1.AND.i1.NE.0.AND.i2.NE.0) THEN
376 C EDGES WHICH DO NOT BELONG TO ANY CELL ARE ALLOWED
377  res=5
378  WRITE(gridman_unit,*)
379  w " GRIDMAN_GRID_CHECK:",
380  w " duplicate cell index for an edge ",
381  w " IEDGE ,CELLS, ",iedge,grid%CELLS(:,iedge)
382  GOTO 100
383  END IF
384  END DO
385 
386 C DOES NCELLS CORRESPOND TO THE MAXIMUM INDEX IN CELLS?
387  ncells_max=0
388  DO iedge=1,grid%NEDGES
389  icell=grid%CELLS(1,iedge)
390  IF(icell.GT.ncells_max) ncells_max=icell
391  icell=grid%CELLS(2,iedge)
392  IF(icell.GT.ncells_max) ncells_max=icell
393  END DO
394  IF(ncells_max.NE.grid%NCELLS) THEN
395  res=6
396  WRITE(gridman_unit,*)
397  w " GRIDMAN_GRID_CHECK: NCELLS and CELLS do not match"
398  WRITE(gridman_unit,*) "NCELLS, NCELLS_MAX ",
399  w grid%NCELLS,ncells_max
400  GOTO 100
401  END IF
402 
403 C CHECK INDICES OF POINTS
404  DO iedge=1,grid%NEDGES
405  DO ie=1,grid%EDIM
406  IF(grid%POINTS(ie,iedge).GT.grid%NPOINTS.OR.
407  f grid%POINTS(ie,iedge).LT.1) THEN
408  res=7 !POINT INDEX IS OUT OF BOUNDS
409  WRITE(gridman_unit,*)
410  w " GRIDMAN_GRID_CHECK: point index is out of bounds ",
411  w " IEDGE, IE, POINTS, NPOINTS ",
412  w iedge, ie, grid%POINTS(ie,iedge),grid%NPOINTS
413  GOTO 100
414  END IF
415  END DO
416  END DO
417 
418 C UNIT SCALE FACTORS
419  IF(grid%UNIT2SI.LT.tiny(grid%UNIT2SI)) THEN
420  res=8
421  WRITE(gridman_unit,*)
422  w " GRIDMAN_GRID_CHECK: incorrect unit scale"
423  WRITE(gridman_unit,*) " GRID%UNIT2SI ",grid%UNIT2SI
424  GOTO 100
425  END IF
426 
427 C EDGE INDICES
428  IF(ALLOCATED(grid%EDGEINDEX)) THEN
429  IF(grid%NEDGEINDEX.LT.1) THEN
430  res=9
431  WRITE(gridman_unit,*) " GRIDMAN_GRID_CHECK: ",
432  w "incorrect number of the groups of edge indices"
433  WRITE(gridman_unit,*) " NEDGEINDEX ",grid%NEDGEINDEX
434  GOTO 100
435  END IF
436  n1=SIZE(grid%EDGEINDEX)
437  IF(n1.NE.grid%NEDGEINDEX) THEN
438  res=9
439  WRITE(gridman_unit,*) " GRIDMAN_GRID_CHECK: ",
440  w "array size mismatch, EDGEINDEX"
441  WRITE(gridman_unit,*) " NEDGEINDEX, SIZE(EDEINDEX) ",
442  w grid%NEDGEINDEX, n1
443  GOTO 100
444  END IF
445  DO id=1,grid%NEDGEINDEX
446  CALL gridman_index_check(grid%EDGEINDEX(id),res,ierr)
447  IF(res.NE.0.OR.ierr.NE.0) THEN
448  res=9
449  WRITE(gridman_unit,*) " GRIDMAN_GRID_CHECK: ",
450  w "incorrect index object ",id
451  GOTO 100
452  END IF
453  DO ie=1,grid%EDGEINDEX(id)%NELEMENTS
454  iedge=grid%EDGEINDEX(id)%INDEXES(0,ie)
455  IF(iedge.LT.0.OR.iedge.GT.grid%NEDGES) THEN
456  res=9
457  WRITE(gridman_unit,*) " GRIDMAN_GRID_CHECK: ",
458  w "incorrect index object ",id
459  WRITE(gridman_unit,*) " Edge index is out of bounds"
460  WRITE(gridman_unit,*) " IE, IEDGE, NEDGES ",
461  w ie, iedge, grid%NEDGES
462  GOTO 100
463  END IF
464  END DO
465  END DO
466  ELSE
467  IF(grid%NEDGEINDEX.GT.1) THEN
468  res=9
469  WRITE(gridman_unit,*) " GRIDMAN_GRID_CHECK: "
470  WRITE(gridman_unit,*) " edge indices are not allocated, ",
471  w " but NEDGEINDEX>1 ",grid%NEDGEINDEX
472  GOTO 100
473  END IF
474  END IF
475 
476 C CELL INDICES
477  IF(ALLOCATED(grid%CELLINDEX)) THEN
478  IF(grid%NCELLINDEX.LT.1) THEN
479  res=10
480  WRITE(gridman_unit,*) " GRIDMAN_GRID_CHECK: ",
481  w "incorrect number of the groups of cell indices"
482  WRITE(gridman_unit,*) " NCELLINDEX ",grid%NCELLINDEX
483  GOTO 100
484  END IF
485  n1=SIZE(grid%CELLINDEX)
486  IF(n1.NE.grid%NCELLINDEX) THEN
487  res=10
488  WRITE(gridman_unit,*) " GRIDMAN_GRID_CHECK: ",
489  w "array size mismatch, CELLINDEX"
490  WRITE(gridman_unit,*) " NCELLINDEX, SIZE(CELLINDEX) ",
491  w grid%NCELLINDEX, n1
492  GOTO 100
493  END IF
494  DO id=1,grid%NCELLINDEX
495  CALL gridman_index_check(grid%CELLINDEX(id),res,ierr)
496  IF(res.NE.0.OR.ierr.NE.0) THEN
497  res=10
498  WRITE(gridman_unit,*) " GRIDMAN_GRID_CHECK: ",
499  w "incorrect index object ",id
500  GOTO 100
501  END IF
502  DO ie=1,grid%CELLINDEX(id)%NELEMENTS
503  icell=grid%CELLINDEX(id)%INDEXES(0,ie)
504  IF(icell.LT.0.OR.icell.GT.grid%NCELLS) THEN
505  res=10
506  WRITE(gridman_unit,*) " GRIDMAN_GRID_CHECK: ",
507  w "incorrect index object ",id
508  WRITE(gridman_unit,*) " Cell index is out of bounds"
509  WRITE(gridman_unit,*) " IE, ICELL, NCELLS ",
510  w ie, icell, grid%NCELLS
511  GOTO 100
512  END IF
513  END DO
514  END DO
515  ELSE
516  IF(grid%NCELLINDEX.GT.1) THEN
517  res=10
518  WRITE(gridman_unit,*) " GRIDMAN_GRID_CHECK: "
519  WRITE(gridman_unit,*) " edge indices are not allocated, ",
520  w "but NCELLINDEX>1 ",grid%NCELLINDEX
521  GOTO 100
522  END IF
523  END IF
524 
525  res=0
526 
527  100 CONTINUE
528 
529  IF(gridman_dbg)
530  w WRITE(gridman_unit,*) "GRIDMAN_GRID_CHECK finished"
531 
532  END SUBROUTINE gridman_grid_check
533 
534 C*****************************************************************************************
535 C> Save grid object in a file defined by the unit index
536 C*****************************************************************************************
537 C>
538 C> GRIDMAN_GRID_UNITWRITE and GRIDMAN_GRID_FILEWRITE
539 C> are combined in the interface GRIDMAN_GRID_WRITE
540  SUBROUTINE gridman_grid_unitwrite(GRID,IOUT,IERR)
542  u gridman_dbg,gridman_check,
543  u gridman_length,gridman_ver
545  IMPLICIT NONE
546  INTRINSIC len_trim
547 C> Grid object to be saved
548  TYPE(gridman_grid) :: GRID
549 C> Index of the output unit
550  INTEGER,INTENT(IN) :: IOUT
551 C> Error code
552  INTEGER,INTENT(OUT) :: IERR
553 
554  CHARACTER*32 :: FRM !FORMAT STRING
555  INTEGER(GRIDMAN_SP) :: IEDGE,IPOINT
556  INTEGER :: IO,II,RES
557 
558  IF(gridman_dbg)
559  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_UNITWRITE"
560 
561  ierr=0
562 
563  IF(gridman_check) THEN
564  CALL gridman_grid_check(grid,res,ierr)
565  IF(res.NE.0.OR.ierr.GT.0) THEN
566  ierr=100
567  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_UNITWRITE: ",
568  w "incorrect grid object"
569  WRITE(gridman_unit,*) " Writing is skipped"
570  RETURN
571  END IF
572  END IF
573 
574 C> At the moment only ASCII format is supported
575 C> @todo HDF5 format has shell be added in the future
576  WRITE(iout,'(A)',iostat=io) 'DATA-TYPE GRIDMAN_GRID'
577  IF(io.NE.0) GOTO 100
578  WRITE(iout,'(A)',iostat=io) 'VERSION'
579  IF(io.NE.0) GOTO 100
580  WRITE(iout,'(I7)',iostat=io) gridman_ver
581  IF(io.NE.0) GOTO 100
582  WRITE(iout,'(A)',iostat=io) "GRIDMAN_LENGTH"
583  IF(io.NE.0) GOTO 100
584  WRITE(iout,'(I6)',iostat=io) gridman_length
585  IF(io.NE.0) GOTO 100
586  WRITE(iout,'(A)',iostat=io) 'DESCRIPTION'
587  IF(io.NE.0) GOTO 100
588  WRITE(frm,210) gridman_length
589  210 FORMAT('(A',i4,')')
590  WRITE(iout,frm,iostat=io) grid%DESCRIPTION
591  IF(io.NE.0) GOTO 100
592 
593  IF(gridman_dbg)
594  f WRITE(gridman_unit,*)
595  w " GRIDMAN_GRID_UNITWRITE. Header finished"
596 
597 C THIS PART IS INCLUDED FOR COMPATIBILITY WITH CODES WHICH DO NOT USE GRID OBJECT
598  WRITE(iout,'(A)',iostat=io) 'EDIM, PDIM'
599  IF(io.NE.0) GOTO 100
600  WRITE(iout,'(2I7)',iostat=io) grid%EDIM,grid%PDIM
601 C
602  WRITE(iout,'(A)',iostat=io) 'TYPE, NEDGES, NCELLS, NPOINTS'
603  IF(io.NE.0) GOTO 100
604  WRITE(iout,'(4I7)',iostat=io)
605  w grid%TYPE,grid%NEDGES,grid%NCELLS,grid%NPOINTS
606  IF(io.NE.0) GOTO 100
607 
608  WRITE(iout,'(A)',iostat=io) 'NEDGEINDEX, NCELLINDEX'
609  IF(io.NE.0) GOTO 100
610  WRITE(iout,'(2I7)',iostat=io) grid%NEDGEINDEX,grid%NCELLINDEX
611  IF(io.NE.0) GOTO 100
612 
613  IF(gridman_dbg)
614  f WRITE(gridman_unit,*)
615  w " GRIDMAN_GRID_UNITWRITE. Dimensions finished"
616 
617  WRITE(iout,'(A)',iostat=io) 'CELLS(1,:), CELLS(2,:)'
618  IF(io.NE.0) GOTO 100
619  DO iedge=1,grid%NEDGES
620  WRITE(iout,'(2I7)',iostat=io) grid%CELLS(1:2,iedge)
621  IF(io.NE.0) GOTO 100
622  END DO
623 
624  WRITE(iout,'(A)',iostat=io) 'POINTS'
625  IF(io.NE.0) GOTO 100
626  WRITE(frm,220) grid%EDIM
627  220 FORMAT('(',i4,'I7)')
628  DO iedge=1,grid%NEDGES
629  WRITE(iout,frm,iostat=io) grid%POINTS(1:grid%EDIM,iedge)
630  IF(io.NE.0) GOTO 100
631  END DO
632 
633  IF(gridman_dbg)
634  f WRITE(gridman_unit,*)
635  w " GRIDMAN_GRID_UNITWRITE. Topology finished"
636 
637  WRITE(iout,'(A)',iostat=io) 'UNITS'
638  IF(io.NE.0) GOTO 100
639  WRITE(frm,230) gridman_length
640  230 FORMAT('(E20.12,1X,A',i4,')')
641  WRITE(iout,frm,iostat=io) grid%UNIT2SI,grid%UNITS
642 
643  WRITE(iout,'(A)',iostat=io) 'X'
644  IF(io.NE.0) GOTO 100
645  WRITE(frm,240) grid%EDIM
646  240 FORMAT('(',i4,'E16.8)')
647  DO ipoint=1,grid%NPOINTS
648  WRITE(iout,frm,iostat=io) grid%X(1:grid%PDIM,ipoint)
649  IF(io.NE.0) GOTO 100
650  END DO
651 
652  IF(gridman_dbg)
653  f WRITE(gridman_unit,*)
654  w " GRIDMAN_GRID_UNITWRITE. Coordinates finished"
655 
656  DO ii=1,grid%NEDGEINDEX
657  IF(.NOT.ALLOCATED(grid%EDGEINDEX(ii)%INDEXES)) THEN
658  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_UNITWRITE: "
659  WRITE(gridman_unit,*) "Group of edge incices is not allocated"
660  WRITE(gridman_unit,*) " II ",ii
661  ierr=100
662  END IF
663  CALL gridman_index_write(grid%EDGEINDEX(ii),iout,ierr)
664  IF(ierr.NE.0) GOTO 1000
665  END DO
666 
667  DO ii=1,grid%NCELLINDEX
668  IF(.NOT.ALLOCATED(grid%CELLINDEX(ii)%INDEXES)) THEN
669  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_UNITWRITE: "
670  WRITE(gridman_unit,*) "Group of cell incices is not allocated"
671  WRITE(gridman_unit,*) " II ",ii
672  ierr=100
673  END IF
674  CALL gridman_index_write(grid%CELLINDEX(ii),iout,ierr)
675  IF(ierr.NE.0) GOTO 1000
676  END DO
677 
678  WRITE(iout,'(A)',iostat=io) "END OF DATA-TYPE GRIDMAN_GRID"
679  IF(io.NE.0) GOTO 100
680 
681  IF(gridman_dbg)
682  w WRITE(gridman_unit,*) "GRIDMAN_GRID_UNITWRITE finished"
683 
684  RETURN
685 
686  100 ierr=300
687  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_UNITWRITE: ",
688  w "cannot write to the file"
689  WRITE(gridman_unit,*) " Writing is skipped"
690  RETURN
691 
692 1000 WRITE(gridman_unit,*) "GRIDMAN_GRID_UNITWRITE terminated"
693 
694  END SUBROUTINE gridman_grid_unitwrite
695 
696 C*****************************************************************************************
697 C> Save grid object in a file defined by the file name
698 C*****************************************************************************************
699 C>
700 C> This subroutine is a shell for GRIDMAN_GRID_UNITWRITE
701 C>
702 C> GRIDMAN_GRID_UNITWRITE and GRIDMAN_GRID_FILEWRITE
703 C> are combined in the interface GRIDMAN_GRID_WRITE
704  SUBROUTINE gridman_grid_filewrite(GRID,FNAME,IERR)
707  IMPLICIT NONE
708  INTRINSIC trim
709 C> Grid object to be saved
710  TYPE(gridman_grid) :: GRID
711 C> Name of the output file
712  CHARACTER(*) :: FNAME
713 C> Error code
714  INTEGER,INTENT(OUT) :: IERR
715  INTEGER :: IO
716  LOGICAL :: LE1,LE2,LO
717 
718  IF(gridman_dbg)
719  w WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_FILEWRITE"
720 
721  ierr=0
722 
723  INQUIRE(file='fort.3',exist=le1)
724  INQUIRE(file='FORT.3',exist=le2)
725  IF(le1.OR.le2) THEN
726  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_FILEWRITE: ",
727  w "can not attach file to UNIT=3 " ,
728  w "because file fort.3 already exists"
729  ierr=300
730  RETURN
731  END IF
732 
733  INQUIRE(unit=3,opened=lo)
734  IF(lo) THEN
735  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_FILEWRITE: ",
736  w "can not attach file to UNIT=3 ",
737  w "because this unit already opened"
738  ierr=300
739  RETURN
740  END IF
741 
742  OPEN(unit=3,file=trim(fname),status='REPLACE',iostat=io)
743  IF(io.NE.0) THEN
744  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_FILEWRITE: ",
745  w "cannot open file ",trim(fname)
746  ierr=300
747  RETURN
748  END IF
749 
750  CALL gridman_grid_unitwrite(grid,3,ierr)
751  CLOSE(3)
752 
753  IF(gridman_dbg)
754  w WRITE(gridman_unit,*) "GRIDMAN_GRID_FILEWRITE finished"
755 
756  END SUBROUTINE gridman_grid_filewrite
757 
758 C*****************************************************************************************
759 C> Read grid object from file defined by unit number
760 C*****************************************************************************************
761 C>
762 C> WARNING: Object which already exists will be overwritten!
763 C>
764 C> GRIDMAN_GRID_UNITREAD and GRIDMAN_GRID_FILEREAD are
765 C> Combined in the interafce GRIDMAN_GRID_READ
766  SUBROUTINE gridman_grid_unitread(GRID,IIN,IERR)
769  u gridman_dbg,gridman_check,gridman_length
772  IMPLICIT NONE
773  INTRINSIC min
774 C> Resulting grid object
775  TYPE(gridman_grid) :: GRID
776 C> Index of the input unit
777  INTEGER,INTENT(IN) :: IIN
778 C> Error code
779  INTEGER,INTENT(OUT) :: IERR
780 
781  CHARACTER*32 :: FRM
782  CHARACTER*128 :: SBUFF
783  CHARACTER(LEN=GRIDMAN_LENGTH) :: DESCRIPTION
784  INTEGER :: IO,LENGTH,EDIM,PDIM,NTYPE,
785  i nedgeindex,ncellindex,ii,ierr0,res
786  INTEGER(GRIDMAN_SP) :: VER,NEDGES,NCELLS,NPOINTS,IEDGE,IPOINT
787 
788  IF(gridman_dbg) WRITE(gridman_unit,*)
789  f "Starting GRIDMAN_GRID_UNITREAD"
790 
791  ierr=0
792 
793 C> At the moment only ASCII format is supported
794 C> @todo HDF5 format has shell be added in the future
795 
796  READ(iin,*,iostat=io) sbuff !'DATA-TYPE GRIDMAN_GRID'
797  IF(io.NE.0) GOTO 100
798  READ(iin,*,iostat=io) sbuff !'VERSION'
799  IF(io.NE.0) GOTO 100
800  READ(iin,'(I7)',iostat=io) ver
801  IF(io.NE.0) GOTO 100
802  READ(iin,*,iostat=io) sbuff !GRIDMAN_LENGTH
803  IF(io.NE.0) GOTO 100
804  READ(iin,'(I6)',iostat=io) length
805  IF(io.NE.0) GOTO 100
806  length=min(length,gridman_length)
807  READ(iin,*,iostat=io) sbuff !'DESCRIPTION'
808  IF(io.NE.0) GOTO 100
809  WRITE(frm,210) length
810  210 FORMAT('(A',i4,')')
811  READ(iin,frm,iostat=io) description
812  IF(io.NE.0) GOTO 100
813 
814  IF(gridman_dbg)
815  f WRITE(gridman_unit,*) " GRIDMAN_GRID_UNITREAD. Header finished"
816 
817  READ(iin,*,iostat=io) sbuff !EDIM,PDIM
818  READ(iin,'(2I7)',iostat=io) edim,pdim
819  IF(io.NE.0) GOTO 100
820 
821  READ(iin,*,iostat=io) sbuff !'TYPE, NEDGES, NCELLS, NPOINTS, LEDGES'
822  IF(io.NE.0) GOTO 100
823  READ(iin,'(4I7)',iostat=io) ntype,nedges,ncells,npoints
824  IF(io.NE.0) GOTO 100
825 
826  READ(iin,*,iostat=io) sbuff !'NEDGEINDEX, NCELLINDEX'
827  IF(io.NE.0) GOTO 100
828  READ(iin,'(2I7)',iostat=io) nedgeindex,ncellindex
829  IF(io.NE.0) GOTO 100
830 
831  IF(gridman_dbg)
832  f WRITE(gridman_unit,*)
833  w " GRIDMAN_GRID_UNITREAD. Dimensions finished"
834 
835  CALL gridman_grid_allocate(grid,ntype,nedges,npoints,ncells,ierr,
836  c nedgeindex,ncellindex)
837  IF(ierr.NE.0) GOTO 1000
838 
839  grid%DESCRIPTION=description
840 
841  READ(iin,*,iostat=io) sbuff !CELLS
842  IF(io.NE.0) GOTO 100
843  DO iedge=1,grid%NEDGES
844  READ(iin,'(3I7)',iostat=io) grid%CELLS(1:2,iedge)
845  IF(io.NE.0) GOTO 100
846  END DO
847 
848  READ(iin,*,iostat=io) sbuff !POINTS
849  IF(io.NE.0) GOTO 100
850  WRITE(frm,220) grid%EDIM
851  220 FORMAT('(',i4,'I7)')
852 
853  DO iedge=1,grid%NEDGES
854  READ(iin,frm,iostat=io) grid%POINTS(1:grid%EDIM,iedge)
855  IF(io.NE.0) GOTO 100
856  END DO
857 
858  IF(gridman_dbg)
859  f WRITE(gridman_unit,*)
860  w " GRIDMAN_GRID_UNITREAD. Topology finished"
861 
862  READ(iin,*,iostat=io) sbuff !UNITS
863  IF(io.NE.0) GOTO 100
864  WRITE(frm,230) length
865  230 FORMAT('(E20.12,1X,A',i4,')')
866  READ(iin,frm,iostat=io) grid%UNIT2SI,grid%UNITS
867  IF(io.NE.0) GOTO 100
868 
869  READ(iin,*,iostat=io) sbuff !X
870  IF(io.NE.0) GOTO 100
871  WRITE(frm,240) grid%PDIM
872  240 FORMAT('(',i4,'E16.8)')
873  DO ipoint=1,grid%NPOINTS
874  READ(iin,frm,iostat=io) grid%X(1:grid%PDIM,ipoint)
875  IF(io.NE.0) GOTO 100
876  END DO
877 
878  IF(gridman_dbg)
879  f WRITE(gridman_unit,*)
880  w " GRIDMAN_GRID_UNITWRITE. Coordinates finished"
881 
882  DO ii=1,grid%NEDGEINDEX
883  CALL gridman_index_read(grid%EDGEINDEX(ii),iin,ierr)
884  IF(ierr.NE.0) GOTO 1000
885  END DO
886 
887  DO ii=1,grid%NCELLINDEX
888  CALL gridman_index_read(grid%CELLINDEX(ii),iin,ierr)
889  IF(ierr.NE.0) GOTO 1000
890  END DO
891 
892  READ(iin,*,iostat=io) sbuff ! END OF DATA-TYPE GRIDMAN_GRID
893  IF(io.NE.0) GOTO 100
894 
895  IF(gridman_check) THEN
896  CALL gridman_grid_check(grid,res,ierr)
897  IF(res.NE.0.OR.ierr.GT.0) THEN
898  ierr=100
899  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_UNITREAD: ",
900  w "incorrect resulting grid"
901  WRITE(gridman_unit,*) " Writing is skipped"
902  CALL gridman_grid_deallocate(grid,ierr0)
903  RETURN
904  END IF
905  END IF
906 
907  IF(gridman_dbg)
908  w WRITE(gridman_unit,*) "GRIDMAN_GRID_UNITREAD finished"
909 
910  RETURN
911 
912  100 ierr=300
913  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_UNITREAD: ",
914  w "can not read data"
915  RETURN
916  1000 WRITE(gridman_unit,*) "GRIDMAN_GRID_UNITREAD is terminated"
917 
918  END SUBROUTINE gridman_grid_unitread
919 
920 C*****************************************************************************************
921 C> Read grid object from file defined by file name
922 C*****************************************************************************************
923 C>
924 C> This subroutine is a shell for GRIDMAN_GRID_UNITREAD
925 C>
926 C> WARNING: Object which already exists will be overwritten!
927 C>
928 C> GRIDMAN_GRID_UNITREAD and GRIDMAN_GRID_FILEREAD are
929 C> Combined in the interafce GRIDMAN_GRID_READ
930  SUBROUTINE gridman_grid_fileread(GRID,FNAME,IERR)
934  IMPLICIT NONE
935  INTRINSIC trim
936 C> Resulting grid object
937  TYPE(gridman_grid) :: GRID
938 C> Name of the output file
939  CHARACTER(*),INTENT(IN) :: FNAME
940 C> Error code
941  INTEGER,INTENT(OUT) :: IERR
942  INTEGER :: IO
943  LOGICAL :: LO
944 
945  IF(gridman_dbg)
946  f WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_FILEREAD"
947 
948  ierr=0
949 
950  INQUIRE(unit=3,opened=lo)
951  IF(lo) THEN
952  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_FILEWRITE: ",
953  w "can not attach file to UNIT=3 ",
954  w "because this unit already opened"
955  ierr=300
956  RETURN
957  END IF
958 
959  OPEN(unit=3,file=trim(fname),status='OLD',iostat=io)
960  IF(io.NE.0) THEN
961  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_FILEREAD: ",
962  w "can not open file ",trim(fname)
963  ierr=300
964  RETURN
965  END IF
966 
967  CALL gridman_grid_unitread(grid,3,ierr)
968  CLOSE(3)
969 
970  IF(gridman_dbg)
971  w WRITE(gridman_unit,*) "GRIDMAN_GRID_FILEREAD finished"
972 
973  END SUBROUTINE gridman_grid_fileread
974 
975 C*****************************************************************************************
976 C> Create a copy of the grid object
977 C*****************************************************************************************
978 C>
979 C> WARNING: If GRID2 already exists it will be overwritten!
980  SUBROUTINE gridman_grid_copy(GRID2,GRID1,IERR)
982  u gridman_dbg,gridman_check
985  IMPLICIT NONE
986 C> Original object
987  TYPE(gridman_grid) :: GRID1
988 C> Copy to be created
989  TYPE(gridman_grid) :: GRID2
990 C> Error code
991  INTEGER,INTENT(OUT) :: IERR
992 
993  INTEGER :: II,RES,IERR0
994 
995  IF(gridman_dbg) WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_COPY"
996 
997  ierr=0
998 
999  IF(gridman_check) THEN
1000  CALL gridman_grid_check(grid1,res,ierr)
1001  IF(res.NE.0.OR.ierr.GT.0) THEN
1002  ierr=100
1003  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_COPY: ",
1004  w "incorrect input grid object"
1005  RETURN
1006  END IF
1007  END IF
1008 
1009  CALL gridman_grid_allocate(grid2,grid1%TYPE,grid1%NEDGES,
1010  s grid1%NPOINTS,grid1%NCELLS,ierr,
1011  s grid1%NEDGEINDEX,grid1%NCELLINDEX)
1012  IF(ierr.NE.0) GOTO 1000
1013 
1014  grid2%DESCRIPTION=grid1%DESCRIPTION
1015  grid2%CELLS=grid1%CELLS
1016  grid2%POINTS=grid1%POINTS
1017  grid2%UNIT2SI=grid1%UNIT2SI
1018  grid2%UNITS=grid1%UNITS
1019  grid2%X=grid1%X
1020 
1021  DO ii=1,grid2%NEDGEINDEX
1022  CALL gridman_index_copy(grid2%EDGEINDEX(ii),
1023  c grid1%EDGEINDEX(ii),ierr)
1024  IF(ierr.NE.0) GOTO 1000
1025  END DO
1026 
1027  DO ii=1,grid2%NCELLINDEX
1028  CALL gridman_index_copy(grid2%CELLINDEX(ii),
1029  c grid1%CELLINDEX(ii),ierr)
1030  IF(ierr.NE.0) GOTO 1000
1031  END DO
1032 
1033  IF(gridman_check) THEN
1034  CALL gridman_grid_check(grid2,res,ierr)
1035  IF(res.NE.0.OR.ierr.GT.0) THEN
1036  ierr=400
1037  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_COPY: ",
1038  w "incorrect resulting grid object"
1039  CALL gridman_grid_deallocate(grid2,ierr0)
1040  RETURN
1041  END IF
1042  END IF
1043 
1044  IF(gridman_dbg) WRITE(gridman_unit,*) "GRIDMAN_GRID_COPY finished"
1045 
1046  RETURN
1047 
1048  1000 WRITE(gridman_unit,*) "GRIDMAN_GRID_COPY is terminated"
1049 
1050  END SUBROUTINE gridman_grid_copy
1051 
1052 C*****************************************************************************************
1053 C> Compare two grid objects
1054 C*****************************************************************************************
1055 C>
1056 C> - RES=0: two grids are equivalent
1057 C> - RES=1: different grid types
1058 C> - RES=2: different dimensions
1059 C> - RES=3: different topology (CELLs)
1060 C> - RES=4: different topology (POINTs)
1061 C> - RES=5: different units
1062 C> - RES=6: different coordinates
1063 C> - RES=7: different description strings
1064 C> - RES=8: different edge indices
1065 C> - RES=9: different cell indices
1066 C
1067  SUBROUTINE gridman_grid_compare(GRID1,GRID2,RES,IERR)
1069  u gridman_check,gridman_dbg
1071  IMPLICIT NONE
1072  INTRINSIC llt,lgt,trim
1073 C> First object to compare
1074  TYPE(gridman_grid) :: GRID1
1075 C> Second object to compare
1076  TYPE(gridman_grid) :: GRID2
1077 C> Code of the result
1078  INTEGER,INTENT(OUT) :: RES
1079 C> Error code
1080  INTEGER,INTENT(OUT) :: IERR
1081 
1082  INTEGER(GRIDMAN_SP) :: IEDGE,IPOINT
1083  INTEGER :: IP,II,RES1
1084 
1085  IF(gridman_dbg)
1086  w WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_COMPARE"
1087 
1088  ierr=0
1089 
1090  IF(gridman_check) THEN
1091  CALL gridman_grid_check(grid1,res,ierr)
1092  IF(res.NE.0.OR.ierr.GT.0) THEN
1093  ierr=100
1094  res=-1
1095  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_COMPARE: ",
1096  w "incorrect GRID1"
1097  RETURN
1098  END IF
1099  CALL gridman_grid_check(grid2,res,ierr)
1100  IF(res.NE.0.OR.ierr.GT.0) THEN
1101  ierr=100
1102  res=-1
1103  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_COMPARE: ",
1104  w "incorrect GRID2"
1105  RETURN
1106  END IF
1107  END IF
1108 
1109 C COMPARING GRID TYPE
1110  IF(grid1%TYPE.NE.grid2%TYPE) THEN
1111  res=1
1112  WRITE(gridman_unit,*) "GRIDMAN_GRID_COMPARE: types are different"
1113  GOTO 100
1114  END IF
1115 
1116 C COMPARING DIMENSIONS
1117  IF((grid1%NEDGES.NE.grid2%NEDGES).OR.
1118  f (grid1%NPOINTS.NE.grid2%NPOINTS).OR.
1119  f (grid1%NCELLS.NE.grid2%NCELLS).OR.
1120  f (grid1%EDIM.NE.grid2%EDIM).OR.
1121  f (grid1%PDIM.NE.grid2%PDIM)) THEN
1122  res=2
1123  WRITE(gridman_unit,*)
1124  w "GRIDMAN_GRID_COMPARE: dimensions are different"
1125  WRITE(gridman_unit,*) " EDIM1, EDIM2 ",
1126  w grid1%EDIM,grid2%EDIM
1127  WRITE(gridman_unit,*) " PDIM1, PDIM2 ",
1128  w grid1%PDIM,grid2%PDIM
1129  WRITE(gridman_unit,*) " NEDGES1, NEDGES2 ",
1130  w grid1%NEDGES,grid2%NEDGES
1131  WRITE(gridman_unit,*) " NPOINTS1, NPOINTS2 ",
1132  w grid1%NPOINTS,grid2%NPOINTS
1133  WRITE(gridman_unit,*) " NCELLS1, NCELLS2 ",
1134  w grid1%NCELLS,grid2%NCELLS
1135  GOTO 100
1136  END IF
1137 
1138 C COMPARING TOPOLOGY
1139  DO iedge=1,grid1%NEDGES
1140  IF( ( grid1%CELLS(1,iedge).NE.grid2%CELLS(1,iedge) ).OR.
1141  f ( grid1%CELLS(2,iedge).NE.grid2%CELLS(2,iedge) ) ) THEN
1142  res=3
1143  WRITE(gridman_unit,*)
1144  w "GRIDMAN_GRID_COMPARE: topology (CELLs) is different"
1145  WRITE(gridman_unit,*) " IEDGE ",iedge
1146  WRITE(gridman_unit,*) " CELL1", grid1%CELLS(:,iedge)
1147  WRITE(gridman_unit,*) " CELL2", grid2%CELLS(:,iedge)
1148  GOTO 100
1149  END IF
1150  DO ip=1,grid1%EDIM
1151  IF(grid1%POINTS(ip,iedge).NE.grid2%POINTS(ip,iedge)) THEN
1152  res=4
1153  WRITE(gridman_unit,*)
1154  w "GRIDMAN_GRID_COMPARE: topology (POINTs) is different"
1155  WRITE(gridman_unit,*) " IEDGE, IP ",iedge,ip
1156  WRITE(gridman_unit,*) " POINT1, POINT2 ",
1157  w grid1%POINTS(ip,iedge),
1158  w grid2%POINTS(ip,iedge)
1159  GOTO 100
1160  END IF
1161  END DO
1162  END DO
1163 
1164 C COMPARING UNITS
1165  IF(diff_real(grid1%UNIT2SI,grid2%UNIT2SI).OR.
1166  f llt(grid1%UNITS,grid2%UNITS).OR.
1167  f lgt(grid1%UNITS,grid2%UNITS) ) THEN
1168  res=5
1169  WRITE(gridman_unit,*)
1170  w "GRIDMAN_GRID_COMPARE: units are different"
1171  WRITE(gridman_unit,*) " UNITS2SI_1, UNIT2SI_2 ",
1172  w grid1%UNIT2SI,grid2%UNIT2SI
1173  WRITE(gridman_unit,*) " UNITS_1, UNITS_2 ",
1174  w trim(grid1%UNITS),trim(grid2%UNITS)
1175  GOTO 100
1176  END IF
1177 
1178 C COMPARING COORDINATES
1179  DO ipoint=1,grid1%NPOINTS
1180  DO ip=1,grid1%PDIM
1181  IF(diff_real(grid1%X(ip,ipoint),grid2%X(ip,ipoint))) THEN
1182  res=6
1183  WRITE(gridman_unit,*)
1184  w "GRIDMAN_GRID_COMPARE: coordinates (X) are different"
1185  WRITE(gridman_unit,*) " IP, IPOINT ",ip,ipoint
1186  WRITE(gridman_unit,*) " X1, X2 ",grid1%X(ip,ipoint),
1187  w grid2%X(ip,ipoint)
1188  GOTO 100
1189  END IF
1190  END DO
1191  END DO
1192 
1193 C COMPARING NAMES (DESCRIPTIONS)
1194  IF(llt(grid1%DESCRIPTION,grid2%DESCRIPTION).OR.
1195  f lgt(grid1%DESCRIPTION,grid2%DESCRIPTION)) THEN
1196  res=7
1197  WRITE(gridman_unit,*)
1198  w "GRIDMAN_GRID_COMPARE: descriptions are different"
1199  WRITE(gridman_unit,*) " DESCRIPTION1: ",trim(grid1%DESCRIPTION)
1200  WRITE(gridman_unit,*) " DESCRIPTION2: ",trim(grid2%DESCRIPTION)
1201  GOTO 100
1202  END IF
1203 
1204 C EDGE INDICES
1205  IF(grid1%NEDGEINDEX.NE.grid2%NEDGEINDEX) THEN
1206  res=8
1207  WRITE(gridman_unit,*) "GRIDMAN_GRID_COMPARE: ",
1208  w "number of the edge index groups is different"
1209  WRITE(gridman_unit,*) " NEDGEINDEX1, NEDGEINDEX2 ",
1210  w grid1%NEDGEINDEX,grid2%NEDGEINDEX
1211  GOTO 100
1212  END IF
1213  DO ii=1,grid1%NEDGEINDEX
1214  CALL gridman_index_compare(grid1%EDGEINDEX(ii),
1215  c grid2%EDGEINDEX(ii),res1,ierr)
1216  IF(res1.NE.0.OR.ierr.NE.0) THEN
1217  WRITE(gridman_unit,*) "GRIDMAN_GRID_COMPARE: ",
1218  w "edge indices are different, II ",ii
1219  res=8
1220  GOTO 100
1221  END IF
1222  END DO
1223 
1224 C CELL INDICES
1225  IF(grid1%NCELLINDEX.NE.grid2%NCELLINDEX) THEN
1226  res=9
1227  WRITE(gridman_unit,*) "GRIDMAN_GRID_COMPARE: ",
1228  w "number of the edge index groups is different"
1229  WRITE(gridman_unit,*) " NCELLINDEX1, NCELLINDEX2 ",
1230  w grid1%NCELLINDEX,grid2%NCELLINDEX
1231  GOTO 100
1232  END IF
1233  DO ii=1,grid1%NCELLINDEX
1234  CALL gridman_index_compare(grid1%CELLINDEX(ii),
1235  c grid2%CELLINDEX(ii),res1,ierr)
1236  IF(res1.NE.0.OR.ierr.NE.0) THEN
1237  WRITE(gridman_unit,*) "GRIDMAN_GRID_COMPARE: ",
1238  w "edge indices are different, II ",ii
1239  res=9
1240  GOTO 100
1241  END IF
1242  END DO
1243 
1244  res=0
1245 
1246  100 CONTINUE
1247  IF(gridman_dbg)
1248  w WRITE(gridman_unit,*) "GRIDMAN_GRID_COMPARE finished"
1249 
1250  CONTAINS
1251 C
1252 C RETURNS .TRUE. IF X1, X2 ARE DIFFERENT AND .FALSE. OTHERWISE
1253 C
1254  FUNCTION diff_real(X1,X2)
1255  USE gridman,ONLY:gridman_dp,gridman_tol
1256  INTRINSIC abs,tiny
1257 
1258  REAL(GRIDMAN_DP),INTENT(IN) :: X1,X2
1259  LOGICAL :: DIFF_REAL
1260  REAL(GRIDMAN_DP) :: EPS
1261 
1262  eps=10.*tiny(x1)
1263  IF(abs(x1-x2).GT.gridman_tol*(abs(x1)+abs(x2)+eps)) THEN
1264  diff_real=.true.
1265  ELSE
1266  diff_real=.false.
1267  END IF
1268 
1269  END FUNCTION diff_real
1270 
1271  END SUBROUTINE gridman_grid_compare
1272 
1273 C*****************************************************************************************
1274 C> Take data from one grid object to another
1275 C*****************************************************************************************
1276 C
1277 C> Opposite to GRIDMAN_GRID_COPY, GRID2 is not overwritten and has to exist.
1278 C> Index tables of GRID2 are not modified
1279 C>
1280 C> WARNING: indices are not taken
1281 C
1282  SUBROUTINE gridman_grid_take(GRID2,GRID1,IERR)
1284  u gridman_dbg,gridman_check
1286  IMPLICIT NONE
1287  INTRINSIC min,ALLOCATED,SIZE
1288 C> Grid object to which the copy is taken
1289  TYPE(gridman_grid) :: GRID2
1290 C> Original grid object
1291  TYPE(gridman_grid) :: GRID1
1292  INTEGER,INTENT(OUT) :: IERR
1293 
1294  INTEGER(GRIDMAN_SP) :: N1,N2
1295  INTEGER :: RES
1296 
1297  IF(gridman_dbg) WRITE(gridman_unit,*) "Starting GRIDMAN_GRID_TAKE"
1298 
1299  IF(gridman_check) THEN
1300  CALL gridman_grid_check(grid1,res,ierr)
1301  IF(res.NE.0.OR.ierr.GT.0) THEN
1302  ierr=100
1303  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_TAKE: ",
1304  w "incorrect grid object"
1305  RETURN
1306  END IF
1307  END IF
1308 
1309  ierr=0
1310 
1311  grid2%DESCRIPTION=grid1%DESCRIPTION
1312  grid2%UNIT2SI=grid1%UNIT2SI
1313  grid2%UNITS=grid1%UNITS
1314 
1315  n1=SIZE(grid2%CELLS,1)
1316  n1=min(2,n1)
1317  n2=SIZE(grid2%CELLS,2)
1318  n2=min(grid1%NEDGES,n2)
1319  IF(ALLOCATED(grid2%CELLS)) THEN
1320  grid2%CELLS(1:n1,1:n2)=grid1%CELLS(1:n1,1:n2)
1321  ELSE
1322  ierr=-100
1323  WRITE(gridman_unit,*) "WARNING from GRIDMAN_GRID_TAKE: ",
1324  w "GRID2%CELLS is not allocated"
1325  END IF
1326 
1327  n1=SIZE(grid2%POINTS,1)
1328  n1=min(grid1%EDIM,n1)
1329  n2=SIZE(grid2%POINTS,2)
1330  n2=min(grid1%NEDGES,n2)
1331  IF(ALLOCATED(grid2%POINTS)) THEN
1332  grid2%POINTS(1:n1,1:n2)=grid1%POINTS(1:n1,1:n2)
1333  ELSE
1334  ierr=-100
1335  WRITE(gridman_unit,*) "WARNING from GRIDMAN_GRID_TAKE: ",
1336  w "GRID2%POINTS is not allocated"
1337  END IF
1338 
1339  n1=SIZE(grid2%X,1)
1340  n1=min(grid1%PDIM,n1)
1341  n2=SIZE(grid2%X,2)
1342  n2=min(grid1%NPOINTS,n2)
1343  IF(ALLOCATED(grid2%X)) THEN
1344  grid2%X(1:n1,1:n2)=grid1%X(1:n1,1:n2)
1345  ELSE
1346  ierr=-100
1347  WRITE(gridman_unit,*) "WARNING from GRIDMAN_GRID_TAKE: ",
1348  w "GRID2%X is not allocated"
1349  END IF
1350 
1351  IF(gridman_dbg) WRITE(gridman_unit,*) "GRIDMAN_GRID_TAKE finished"
1352 
1353  END SUBROUTINE gridman_grid_take
1354 
1355 
1356 C*****************************************************************************************
1357 C> Return the number of cells plus the number of edges not belonging to any cell
1358 C*****************************************************************************************
1359  SUBROUTINE gridman_grid_count(GRID,NCOUNT,IERR)
1361  u gridman_dbg,gridman_check
1363  IMPLICIT NONE
1364  INTRINSIC max
1365 C> Grid object where cells and surfaces are counted
1366  TYPE(gridman_grid) :: GRID
1367 C> Returned number of cells and "free" edges
1368  INTEGER(GRIDMAN_SP),INTENT(OUT) :: NCOUNT
1369 C> Error code
1370  INTEGER,INTENT(OUT) :: IERR
1371 
1372  INTEGER(GRIDMAN_SP) :: IEDGE
1373  INTEGER :: RES
1374 
1375  IF(gridman_dbg) WRITE(gridman_unit,*)
1376  w "Starting GRIDMAN_GRID_COUNT"
1377 
1378  ierr=0
1379 
1380  IF(gridman_check) THEN
1381  CALL gridman_grid_check(grid,res,ierr)
1382  IF(res.NE.0.OR.ierr.GT.0) THEN
1383  ierr=100
1384  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_COUNT: ",
1385  w "incorrect grid object"
1386  RETURN
1387  END IF
1388  END IF
1389 
1390  ncount=max(grid%NCELLS,0)
1391  DO iedge=1,grid%NEDGES
1392  IF(grid%CELLS(1,iedge).LT.1.AND.
1393  f grid%CELLS(2,iedge).LT.1) ncount=ncount+1
1394  END DO
1395 
1396  IF(gridman_dbg) WRITE(gridman_unit,*)
1397  w "GRIDMAN_GRID_COUNT finished"
1398 
1399  END SUBROUTINE gridman_grid_count
1400 
1401 C*****************************************************************************************
1402 C> Print metadata into text array
1403 C>
1404 C> Called from the interface GRIDMAN_GRID_METADATA
1405 C*****************************************************************************************
1406  SUBROUTINE gridman_grid_metadata_text(GRID,TEXT,NT,IERR)
1408  u gridman_dbg,gridman_unit
1409  INTRINSIC sum,trim,minval,maxval
1410 C> Grid whose metadata has to be printed
1411  TYPE(gridman_grid) :: GRID
1412 C> Array of strings with printed metadata, TEXT(NT)
1413  CHARACTER(GRIDMAN_LENGTH),ALLOCATABLE :: TEXT(:)
1414 C> Length of array TEXT
1415  INTEGER,INTENT(OUT) :: NT
1416 C> Error code
1417  INTEGER,INTENT(OUT) :: IERR
1418 
1419  INTEGER :: IS,I,K,J,N
1420 
1421  IF(gridman_dbg) WRITE(gridman_unit,*)
1422  w "Starting GRIDMAN_GRID_METADATA_TEXT"
1423 
1424  ierr=0
1425 
1426  nt=8
1427  IF(grid%NCELLINDEX.GT.0)
1428  f nt=nt+grid%NCELLINDEX*5+sum(grid%CELLINDEX(:)%NINDEX)
1429  IF(grid%NEDGEINDEX.GT.0)
1430  f nt=nt+grid%NEDGEINDEX*5+sum(grid%EDGEINDEX(:)%NINDEX)
1431 
1432  ALLOCATE(text(nt),stat=is)
1433  IF(is.NE.0) THEN
1434  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_METADATA_TEXT: ",
1435  w "cannot allocate array of strings"
1436  WRITE(gridman_unit,*) " NT, LENGTH ",nt,gridman_length
1437  ierr=200
1438  RETURN
1439  END IF
1440 
1441  WRITE(text(1),'(A,A)') 'Description: ',trim(grid%DESCRIPTION)
1442  WRITE(text(2),'(A,I0)') 'Dimensionality: ',grid%PDIM
1443  WRITE(text(3),'(A,I0)') 'Number of cells: ',grid%NCELLS
1444  WRITE(text(4),'(A,I0)') 'Number of edges: ',grid%NEDGES
1445  WRITE(text(5),'(A,I0)') 'Number of points: ',grid%NPOINTS
1446  WRITE(text(6),'(A,A)') 'Units: ',trim(grid%UNITS)
1447  WRITE(text(7),*) 'Minimal values of coordinates: ',
1448  w (minval(grid%X(j,:)),j=1,grid%PDIM)
1449  WRITE(text(8),*) 'Maximal values of coordinates: ',
1450  w (maxval(grid%X(j,:)),j=1,grid%PDIM)
1451 
1452  k=8
1453  DO i=1,grid%NCELLINDEX
1454  k=k+1
1455  WRITE(text(k),'(A,I0)') 'Cell index ',i
1456  k=k+1
1457  WRITE(text(k),'(A,A)') ' Description: ',
1458  w trim(grid%CELLINDEX(i)%DESCRIPTION)
1459  n=grid%CELLINDEX(i)%NINDEX
1460  DO j=1,n
1461  k=k+1
1462  WRITE(text(k),'(A,I0,A,A)') ' Column ',i,': ',
1463  w trim(grid%CELLINDEX(i)%COLUMNS(j))
1464  END DO
1465  k=k+1
1466  WRITE(text(k),'(A,I0)') ' Number of elements: ',
1467  w grid%CELLINDEX(i)%NELEMENTS
1468  k=k+1
1469  WRITE(text(k),*) ' Minimal value(s): ',
1470  w (minval(grid%CELLINDEX(i)%INDEXES(j,:)),j=1,n)
1471  k=k+1
1472  WRITE(text(k),*) ' Maximal value(s): ',
1473  w (maxval(grid%CELLINDEX(i)%INDEXES(j,:)),j=1,n)
1474  END DO
1475 
1476  DO i=1,grid%NEDGEINDEX
1477  k=k+1
1478  WRITE(text(k),'(A,I0)') 'Edge index ',i
1479  k=k+1
1480  WRITE(text(k),'(A,A)') ' Description: ',
1481  w trim(grid%EDGEINDEX(i)%DESCRIPTION)
1482  n=grid%EDGEINDEX(i)%NINDEX
1483  DO j=1,n
1484  k=k+1
1485  WRITE(text(k),'(A,I0,A,A)') ' Column ',i,': ',
1486  w trim(grid%EDGEINDEX(i)%COLUMNS(j))
1487  END DO
1488  k=k+1
1489  WRITE(text(k),'(A,I0)') ' Number of elements: ',
1490  w grid%EDGEINDEX(i)%NELEMENTS
1491  k=k+1
1492  WRITE(text(k),*) ' Minimal value(s): ',
1493  w (minval(grid%EDGEINDEX(i)%INDEXES(j,:)),j=1,n)
1494  k=k+1
1495  WRITE(text(k),*) ' Maximal value(s): ',
1496  w (maxval(grid%EDGEINDEX(i)%INDEXES(j,:)),j=1,n)
1497  END DO
1498 
1499  IF(gridman_dbg) WRITE(gridman_unit,*)
1500  w "GRIDMAN_GRID_METADATA_TEXT finished"
1501 
1502  END SUBROUTINE gridman_grid_metadata_text
1503 
1504 C*****************************************************************************************
1505 C> Print metadata into unit - shell for \ref GRIDMAN_GRID_METADATA_TEXT
1506 C>
1507 C> Called from the interface GRIDMAN_GRID_METADATA
1508 C*****************************************************************************************
1509  SUBROUTINE gridman_grid_metadata_unit(GRID,IUNIT,IERR)
1511  u gridman_dbg,gridman_unit
1513  INTRINSIC trim
1514 C> Grid whose metadata has to be printed
1515  TYPE(gridman_grid) :: GRID
1516 C> Number of the output unit
1517  INTEGER,INTENT(IN) :: IUNIT
1518 C> Error code
1519  INTEGER,INTENT(OUT) :: IERR
1520 
1521  CHARACTER(GRIDMAN_LENGTH),ALLOCATABLE :: TEXT(:)
1522  INTEGER :: I,NT,IS
1523 
1524  IF(gridman_dbg) WRITE(gridman_unit,*)
1525  w "Starting GRIDMAN_GRID_METADATA_UNIT"
1526 
1527  ierr=0
1528 
1529  CALL gridman_grid_metadata_text(grid,text,nt,ierr)
1530  IF(ierr.NE.0) THEN
1531  WRITE(gridman_unit,*) "GRIDMAN_GRID_METADATA_UNIT terminated"
1532  RETURN
1533  END IF
1534 
1535  DO i=1,nt
1536  WRITE(iunit,*,iostat=is) trim(text(i))
1537  IF(is.NE.0) THEN
1538  WRITE(gridman_unit,*) "ERROR in GRIDMAN_GRID_METADATA_UNIT: ",
1539  w "printing failed, line ",i
1540  DEALLOCATE(text)
1541  ierr=300
1542  RETURN
1543  END IF
1544  END DO
1545 
1546  DEALLOCATE(text,stat=is)
1547  IF(is.NE.0) THEN
1548  ierr=-200
1549  WRITE(gridman_unit,*)
1550  w "WARNING from GRIDMAN_GRID_METADATA_UNIT: ",
1551  w "memory deallocation error"
1552  END IF
1553 
1554  IF(gridman_dbg) WRITE(gridman_unit,*)
1555  w "GRIDMAN_GRID_METADATA_UNIT finished"
1556 
1557  END SUBROUTINE gridman_grid_metadata_unit
subroutine gridman_grid_check(GRID, RES, IERR)
Check consistency of the grid data.
Definition: grid1.f:289
subroutine gridman_index_check(INDEX, RES, IERR)
Check index object.
Definition: index.f:149
integer, parameter, public gridman_length
Length of the description strings.
Definition: gridman.f:101
integer, save, public gridman_unit
Index of the standard output unit.
Definition: gridman.f:116
real(gridman_dp), save, public gridman_tol
Tolerance parameter which is used to compare two real numbers.
Definition: gridman.f:127
integer, parameter, public gridman_dp
Kind parameter for real numbers.
Definition: gridman.f:93
Explicit interfaces to GRIDMAN subroutines and functions.
Definition: gridman.f:251
subroutine gridman_grid_allocate(GRID, TYPE, NEDGES, NPOINTS, NCELLS, IERR, NEDGEINDEX, NCELLINDEX)
Allocate GRIDMAN_GRID object.
Definition: grid1.f:30
subroutine gridman_index_compare(INDEX1, INDEX2, RES, IERR)
Compare two index objects.
Definition: index.f:423
subroutine gridman_index_copy(INDEX2, INDEX1, IERR)
Create a copy of the index object.
Definition: index.f:538
subroutine gridman_grid_unitread(GRID, IIN, IERR)
Read grid object from file defined by unit number.
Definition: grid1.f:767
subroutine gridman_grid_take(GRID2, GRID1, IERR)
Take data from one grid object to another.
Definition: grid1.f:1283
subroutine gridman_grid_count(GRID, NCOUNT, IERR)
Return the number of cells plus the number of edges not belonging to any cell.
Definition: grid1.f:1360
Data-type which describes a grid as a set of edges, methods in grid.f.
Definition: gridman.f:168
subroutine gridman_grid_metadata_unit(GRID, IUNIT, IERR)
Print metadata into unit - shell for GRIDMAN_GRID_METADATA_TEXT.
Definition: grid1.f:1510
subroutine gridman_grid_deallocate(GRID, IERR)
Deallocate grid object.
Definition: grid1.f:184
subroutine gridman_grid_copy(GRID2, GRID1, IERR)
Create a copy of the grid object.
Definition: grid1.f:981
subroutine gridman_grid_unitwrite(GRID, IOUT, IERR)
Save grid object in a file defined by the unit index.
Definition: grid1.f:541
logical, save, public gridman_dbg
Switch for debugging mode.
Definition: gridman.f:122
subroutine gridman_grid_metadata_text(GRID, TEXT, NT, IERR)
Print metadata into text array.
Definition: grid1.f:1407
subroutine gridman_grid_fileread(GRID, FNAME, IERR)
Read grid object from file defined by file name.
Definition: grid1.f:931
subroutine gridman_grid_filewrite(GRID, FNAME, IERR)
Save grid object in a file defined by the file name.
Definition: grid1.f:705
subroutine gridman_index_read(INDEX, IIN, IERR)
Read index object from file.
Definition: index.f:327
subroutine gridman_index_write(INDEX, IOUT, IERR)
Save index object in a file.
Definition: index.f:234
Definition of data types, global constants and variables.
Definition: gridman.f:83
subroutine gridman_index_deallocate(INDEX, IERR)
Allocate index object.
Definition: index.f:98
subroutine gridman_grid_compare(GRID1, GRID2, RES, IERR)
Compare two grid objects.
Definition: grid1.f:1068
integer, parameter, public gridman_sp
Kind parameter for integer numbers.
Definition: gridman.f:95