restart.f90
Go to the documentation of this file.
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: restart.f90 #
5 !# #
6 !# Copyright (C) 2015-2019 #
7 !# Manuel Jung <mjung@astrophysik.uni-kiel.de> #
8 !# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
9 !# #
10 !# This program 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 2 of the License, or (at #
13 !# your option) any later version. #
14 !# #
15 !# This program is distributed in the hope that it will be useful, but #
16 !# WITHOUT ANY WARRANTY; without even the implied warranty of #
17 !# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
18 !# NON INFRINGEMENT. See the GNU General Public License for more #
19 !# details. #
20 !# #
21 !# You should have received a copy of the GNU General Public License #
22 !# along with this program; if not, write to the Free Software #
23 !# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
24 !# #
25 !#############################################################################
26 
27 !----------------------------------------------------------------------------!
48 !----------------------------------------------------------------------------!
49 PROGRAM restart
50  USE fosite_mod
51 #ifdef PARALLEL
52 #ifdef HAVE_MPI_MOD
53  USE mpi
54 #endif
55 #endif
56  IMPLICIT NONE
57 #ifdef PARALLEL
58 #ifdef HAVE_MPIF_H
59  include 'mpif.h'
60 #endif
61 #endif
62 #ifdef PARALLEL
63 #define OFFSET_TYPE INTEGER(KIND=MPI_OFFSET_KIND)
64 #else
65 #define OFFSET_TYPE INTEGER
66 #endif
67  !--------------------------------------------------------------------------!
68  CLASS(fosite), ALLOCATABLE :: Sim
69  TYPE(Dict_TYP), POINTER :: input
70  REAL :: time,stoptime
71  INTEGER :: i, step
72  CHARACTER(LEN=100) :: filename, filename_tmp
73 #ifdef PARALLEL
74  INTEGER, DIMENSION(3) :: decomposition
75 #endif
76  LOGICAL :: file_exists
77  REAL :: HRATIO = 0.05
78  REAL, DIMENSION(:,:,:), POINTER :: r, Sigma
79  REAL, DIMENSION(:,:,:,:), POINTER :: r_faces
80  REAL, DIMENSION(:,:,:), POINTER :: bccsound
81  REAL, DIMENSION(:,:,:,:), POINTER :: fcsound
82  REAL, PARAMETER :: MBH1 = 0.4465*1.0
83  REAL, PARAMETER :: MBH2 = 0.4465*1.0
84  INTEGER :: count
85  !--------------------------------------------------------------------------!
86 
87  ! load file
88  CALL get_command_argument(1, filename)
89  INQUIRE(file=trim(filename), exist=file_exists)
90  IF(.NOT.file_exists) THEN
91  print *, "restart, ", "Input file does not exist!"
92  stop
93  END IF
94 
95  ALLOCATE(sim)
96 
97  CALL sim%InitFosite()
98 
99  !-------------------- load configuration-----------------------------------!
100  input => loadconfig(trim(filename))
101  CALL getattr(input, "config", sim%config)
102 
103  ! set starting time in fosite
104  CALL getattr(input, '/timedisc/time', time)
105 !time = 0.0
106  CALL setattr(sim%config,"/timedisc/starttime", time)
107 
108 #ifdef PARALLEL
109  decomposition(1) = 1
110  decomposition(2) = -1
111  decomposition(3) = 1
112  CALL setattr(sim%config, "/mesh/decomposition", decomposition)
113 #endif
114 
115  ! get step number by stripping from filename
116  filename_tmp = filename
117  i = scan(filename_tmp,'_',back=.true.)
118  IF(i.GT.0) THEN
119  filename_tmp = filename_tmp(i+1:i+4)
120  IF(filename_tmp(1:1) .EQ. "0") filename_tmp = filename_tmp(2:)
121  IF(filename_tmp(1:1) .EQ. "0") filename_tmp = filename_tmp(2:)
122  IF(filename_tmp(1:1) .EQ. "0") filename_tmp = filename_tmp(2:)
123  READ(filename_tmp,*) step
124  CALL setattr(sim%config,"/datafile/step", step)
125  ELSE
126  CALL sim%Error("Restart", "Unable to find starting step in filename_tmp.")
127  END IF
128 
129 
130  !-------------------- setup & data initialization -------------------------!
131  CALL sim%Setup()
132 
133  CALL loaddata(sim,trim(filename))
134 
135  !--------------------------------------------------------------------------!
136  SELECT TYPE (phys => sim%Physics)
137  TYPE IS(physics_eulerisotherm)
138  ! Set sound speed
139  ALLOCATE( &
140  bccsound(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX,sim%Mesh%KGMIN:sim%Mesh%KGMAX), &
141  fcsound(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX,sim%Mesh%KGMIN:sim%Mesh%KGMAX,6))
142 
143  r => sim%Mesh%RemapBounds(sim%Mesh%radius%bcenter)
144  r_faces => sim%Mesh%RemapBounds(sim%Mesh%radius%faces)
145  bccsound = hratio*sqrt((mbh1+mbh2)*phys%Constants%GN/r(:,:,:))
146  DO i =1,6
147  fcsound(:,:,:,i) = hratio*sqrt((mbh1+mbh2)*phys%constants%GN/r_faces(:,:,:,i))
148  END DO
149  ! set isothermal sound speeds
150  CALL phys%SetSoundSpeeds(sim%Mesh,bccsound)
151  CALL phys%SetSoundSpeeds(sim%Mesh,fcsound)
152 
153  ! boundary conditions
154  ! custom boundary conditions at western boundary if requested
155  SELECT TYPE(bwest => sim%Timedisc%Boundary%boundary(west)%p)
156  CLASS IS (boundary_custom)
157  CALL bwest%SetCustomBoundaries(sim%Mesh,sim%Physics, &
158  (/custom_nograd,custom_outflow,custom_kepler/))
159  END SELECT
160  SELECT TYPE(beast => sim%Timedisc%Boundary%boundary(east)%p)
161  CLASS IS (boundary_custom)
162  CALL beast%SetCustomBoundaries(sim%Mesh,sim%Physics, &
163  (/custom_nograd,custom_outflow,custom_kepler/))
164  END SELECT
165  END SELECT
166  !--------------------------------------------------------------------------!
167 
168 
169  CALL sim%Info(" DATA-----> initial condition: restarted from data file - " // &
170  trim(filename))
171 
172  CALL sim%Physics%Convert2Conservative(sim%Timedisc%pvar,sim%Timedisc%cvar)
173 ! CALL Sim%FirstStep()
174 
175  CALL sim%Run()
176  CALL sim%Finalize()
177  DEALLOCATE(sim)
178 
179 CONTAINS
180 
181 
186 FUNCTION loadconfig(filename) RESULT(res)
187  IMPLICIT NONE
188  !--------------------------------------------------------------------------!
189  CHARACTER(LEN=*) :: filename
190  TYPE(Dict_TYP),POINTER :: res
191  !--------------------------------------------------------------------------!
192  INTEGER :: file, error, keylen, intsize, realsize, &
193  type, bytes, l, dims(5)
194  CHARACTER(LEN=6) :: magic
195  CHARACTER(LEN=2) :: endian
196  CHARACTER(LEN=13) :: header
197  CHARACTER(LEN=4) :: sizestr
198  CHARACTER(LEN=1),DIMENSION(128) :: buffer
199 #ifndef PARALLEL
200  INTEGER :: offset
201  CHARACTER(LEN=1) :: version
202 #else
203  INTEGER(KIND=MPI_OFFSET_KIND) :: offset, offset_0, filesize
204  INTEGER :: version
205 #endif
206  CHARACTER(LEN=1),DIMENSION(:),POINTER :: keybuf
207  CHARACTER(LEN=MAX_CHAR_LEN) :: key,kf
208  REAL,DIMENSION(:,:,:), POINTER :: ptr3 => null()
209  CHARACTER(LEN=1),DIMENSION(:),POINTER :: val
210 #ifdef PARALLEL
211  INTEGER :: handle
212  INTEGER :: ierror
213  INTEGER :: bufsize
214  INTEGER :: position
215  INTEGER, DIMENSION(2) :: gsizes,lsizes,indices,memsizes
216  INTEGER, DIMENSION(MPI_STATUS_SIZE) :: status
217 #endif
218  !--------------------------------------------------------------------------!
219  INTENT(IN) :: filename
220  !--------------------------------------------------------------------------!
221  NULLIFY(res)
222  offset = 1
223 #ifndef PARALLEL
224  file = 5555
225  OPEN(file, &
226  file = trim(filename), &
227  status = 'OLD', &
228  access = 'STREAM', &
229  action = 'READ', &
230  position = 'REWIND', &
231  iostat = error)
232 #else
233  ! Open File
234  CALL mpi_file_open(mpi_comm_world,trim(filename),mpi_mode_rdonly, &
235  mpi_info_null,handle,error)
236  offset_0 = 0
237  CALL mpi_file_set_view(handle,offset_0,mpi_byte,mpi_byte,'native',mpi_info_null,ierror)
238  CALL mpi_file_seek(handle,offset-1,mpi_seek_set,ierror)
239 #endif
240 
241 #ifndef PARALLEL
242  READ(file,pos=offset) magic, endian, version, sizestr
243 #else
244  CALL mpi_file_read_all(handle, header, len(header), mpi_byte, &
245  status,ierror)
246  WRITE (magic, '(A6)') header(1:6)
247  WRITE (endian, '(A2)') header(7:8)
248  version = iachar(header(9:9))
249  WRITE (sizestr, '(A4)') header(10:13)
250 #endif
251  offset = offset + 13
252  READ(sizestr, '(I2,I2)') realsize, intsize
253 
254 #ifndef PARALLEL
255  READ(file, iostat=error, pos=offset) keylen
256 #else
257 ! TODO: keylen from data does not need to have the same intsize like from the actual run with Fosite
258  buffer = ''
259  CALL mpi_file_seek(handle,offset-1,mpi_seek_set,ierror)
260  CALL mpi_file_read_all(handle, buffer, intsize, mpi_byte, status,ierror)
261  keylen = transfer(buffer(1:intsize),keylen)
262 #endif
263  offset = offset + intsize
264 
265  DO WHILE(error.EQ.0)
266  key = ''
267  buffer = ''
268  ALLOCATE(keybuf(keylen))
269 #ifndef PARALLEL
270  READ(file, pos=offset) keybuf,type,bytes
271 #else
272  CALL mpi_file_seek(handle,offset-1,mpi_seek_set,ierror)
273  CALL mpi_file_read_all(handle, buffer, keylen+2*intsize, mpi_byte, &
274  status,ierror)
275  keybuf = transfer(buffer(1:keylen), keybuf)
276  type = transfer(buffer(keylen+1:keylen+intsize), type)
277  bytes = transfer(buffer(keylen+intsize+1:keylen+2*intsize), bytes)
278 #endif
279  WRITE(kf,'(A, I4, A)') '(',keylen,'(A))'
280  WRITE(key,kf) keybuf
281  DEALLOCATE(keybuf)
282  key = trim(key)
283  offset = offset + keylen + 2*intsize
284  dims(:) = 1
285  SELECT CASE(type)
286  CASE(dict_real_twod)
287  l = 3
288  CASE(dict_real_threed)
289  l = 3
290  CASE(dict_real_fourd)
291  l = 4
292  CASE(dict_real_fived)
293  l = 5
294  CASE DEFAULT
295  l = 0
296  END SELECT
297  IF(l.GE.2) THEN
298 #ifndef PARALLEL
299  READ(file, pos=offset) dims(1:l)
300 #else
301  buffer = ''
302  CALL mpi_file_seek(handle,offset-1,mpi_seek_set,ierror)
303  CALL mpi_file_read_all(handle,buffer(1:l*intsize),l*intsize,mpi_byte,status,error)
304  dims(1:l) = transfer(buffer(1:l*intsize),dims(1:l))
305 #endif
306  bytes = bytes - l*intsize
307  offset = offset + l*intsize
308  END IF
309 
310 
311  ! Here, the data fields are skipped (only single values are set, e.g.
312  ! configs, central mass, etc.)
313  ! TODO: The allocation of the field is not good here and just a workaround
314  ! - Problem: The data parts in the file need to be skipped somewhow,
315  ! but the POS argument is not available in F95, which can be used
316  ! on NEC/SX-Ace. At the field decomposition by MPI is done after
317  ! reading in the dictionary in SetupFosite.
318  SELECT CASE(l)
319  CASE(2)
320  CASE(3)
321  CASE(4)
322  CASE(5)
323  CASE DEFAULT
324  IF(bytes.GT.0) THEN
325  ALLOCATE(val(bytes))
326 #ifndef PARALLEL
327  READ(file, pos=offset) val
328 #else
329  buffer = ''
330  CALL mpi_file_seek(handle,offset-1,mpi_seek_set,ierror)
331  CALL mpi_file_read_all(handle,buffer(1:bytes),bytes,mpi_byte,status,error)
332  val = transfer(buffer(1:bytes),val)
333 #endif
334  IF(key.EQ.'/config/mesh/decomposition')THEN
335  ! Do not set the key for composition. Use new one.
336  ELSE
337  CALL setattr(res,key,val,type)
338  END IF
339  DEALLOCATE(val)
340  END IF
341  END SELECT
342 
343  offset = offset + bytes
344 
345 #ifndef PARALLEL
346  READ(file, pos=offset, iostat=error) keylen
347 #else
348  buffer = ''
349  CALL mpi_file_seek(handle,offset-1,mpi_seek_set,ierror)
350  CALL mpi_file_read_all(handle,buffer(1:intsize),intsize,mpi_byte,status,ierror)
351  keylen = transfer(buffer(1:intsize),keylen)
352 
353 
354  CALL mpi_file_get_position(handle,offset_0,ierror)
355  CALL mpi_file_get_size(handle,filesize,ierror)
356  IF (filesize.LE.offset_0) THEN
357  ierror=1
358  END IF
359 
360  error = ierror
361 #endif
362  offset = offset + intsize
363 
364  END DO
365 
366  CLOSE(file)
367 END FUNCTION loadconfig
368 
375 SUBROUTINE loaddata(this,filename)
376  IMPLICIT NONE
377  !--------------------------------------------------------------------------!
378  CLASS(fosite) :: this
379  CHARACTER(LEN=*) :: filename
380  !--------------------------------------------------------------------------!
381  INTEGER :: unit, error, keylen, &
382  intsize, realsize, type, bytes, &
383  l, dims(5)
384 ! CHARACTER(LEN=64) :: keybufsize
385  CHARACTER(LEN=13) :: header
386  CHARACTER(LEN=6) :: magic
387  CHARACTER(LEN=2) :: endian
388  CHARACTER(LEN=1),DIMENSION(128) :: buffer
389 #ifndef PARALLEL
390  INTEGER :: offset
391  CHARACTER(LEN=1) :: version
392 #else
393  INTEGER(KIND=MPI_OFFSET_KIND) :: offset, offset_0, filesize
394  INTEGER :: version
395 #endif
396  CHARACTER(LEN=4) :: sizestr
397  CHARACTER(LEN=1),DIMENSION(:),POINTER :: keybuf
398  CHARACTER(LEN=MAX_CHAR_LEN) :: key,kf
399  REAL,DIMENSION(:,:), POINTER :: ptr2 => null()
400  REAL,DIMENSION(:,:,:), POINTER :: ptr3 => null()
401  REAL,DIMENSION(:,:,:,:), POINTER :: ptr4 => null()
402  REAL,DIMENSION(:,:,:,:,:), POINTER :: ptr5 => null()
403  CHARACTER(LEN=1),DIMENSION(:),POINTER :: val
404 #ifdef PARALLEL
405  INTEGER :: handle
406  INTEGER :: filetype
407  INTEGER :: ierror
408  INTEGER :: bufsize
409  INTEGER :: position
410  INTEGER, DIMENSION(2) :: gsizes,lsizes,indices,memsizes
411  INTEGER, DIMENSION(MPI_STATUS_SIZE) :: status
412 #endif
413  CLASS(sources_base), POINTER :: srcptr
414  CLASS(gravity_base), POINTER :: gravptr
415  !--------------------------------------------------------------------------!
416  INTENT(IN) :: filename
417  INTENT(INOUT) :: this
418  !--------------------------------------------------------------------------!
419  offset = 1
420 
421 #ifndef PARALLEL
422  unit = 5555
423  OPEN(unit, &
424  file = trim(filename), &
425  status = 'OLD', &
426  access = 'STREAM', &
427  action = 'READ', &
428  position = 'REWIND', &
429  iostat = error)
430 #else
431  ! First create Datatype that will be read in
432  gsizes(1) = sim%Mesh%INUM
433  gsizes(2) = sim%Mesh%JNUM
434  lsizes(1) = sim%Mesh%IMAX-sim%Mesh%IMIN+1
435  lsizes(2) = sim%Mesh%JMAX-sim%Mesh%JMIN+1
436  indices(1)= sim%Mesh%IMIN-1
437  indices(2)= sim%Mesh%JMIN-1
438  bufsize = product(lsizes)
439  CALL mpi_type_create_subarray(2, gsizes, lsizes, indices, mpi_order_fortran,&
440  default_mpi_real,filetype,ierror)
441  CALL mpi_type_commit(filetype,ierror)
442 
443  ! Open File
444  CALL mpi_file_open(mpi_comm_world,trim(filename),mpi_mode_rdonly, &
445  mpi_info_null,handle,error)
446 ! CALL MPI_File_seek(handle,offset-1,MPI_SEEK_SET,ierror)
447 #endif
448 
449 #ifndef PARALLEL
450  READ(unit) magic, endian, version, sizestr
451 #else
452  CALL mpi_file_read_all(handle, header, len(header), mpi_byte, &
453  status,ierror)
454  WRITE (magic, '(A6)') header(1:6)
455  WRITE (endian, '(A2)') header(7:8)
456  version = iachar(header(9:9))
457  WRITE (sizestr, '(A4)') header(10:13)
458 #endif
459  offset = offset + 13
460  READ(sizestr, '(I2,I2)') realsize, intsize
461 #ifndef PARALLEL
462  READ(unit, iostat=error) keylen
463 #else
464 ! TODO: keylen from data does not need to have the same intsize like from the actual run with Fosite
465  buffer = ''
466  CALL mpi_file_read_all(handle, buffer, intsize, mpi_byte, status,ierror)
467  keylen = transfer(buffer(1:intsize),keylen)
468 #endif
469  offset = offset + intsize
470 
471 !------------------- loop over data ------------------------------------------!
472  DO WHILE(error.EQ.0)
473  NULLIFY(ptr2)
474  key = ''
475  buffer = ''
476  ALLOCATE(keybuf(keylen))
477 #ifndef PARALLEL
478  READ(unit) keybuf,type,bytes
479 #else
480  CALL mpi_file_read_all(handle, buffer, keylen+2*intsize, mpi_byte, &
481  status,ierror)
482  keybuf = transfer(buffer(1:keylen), keybuf)
483  type = transfer(buffer(keylen+1:keylen+intsize), type)
484  bytes = transfer(buffer(keylen+intsize+1:keylen+2*intsize), bytes)
485 #endif
486  WRITE(kf,'(A, I4, A)') '(',keylen,'(A))'
487  WRITE(key,kf) keybuf
488  DEALLOCATE(keybuf)
489  offset = offset + keylen + 2*intsize
490  dims(:) = 1
491  SELECT CASE(type)
492  CASE(dict_real_twod)
493  l = 3
494  CASE(dict_real_threed)
495  l = 3
496  CASE(dict_real_fourd)
497  l = 4
498  CASE(dict_real_fived)
499  l = 5
500  CASE DEFAULT
501  l = 0
502  END SELECT
503  IF(l.GE.2) THEN
504 #ifndef PARALLEL
505  READ(unit) dims(1:l)
506 #else
507  buffer = ''
508  CALL mpi_file_read_all(handle,buffer(1:l*intsize),l*intsize,mpi_byte,status,error)
509  dims(1:l) = transfer(buffer(1:l*intsize),dims(1:l))
510 #endif
511  bytes = bytes - l*intsize
512  offset = offset + l*intsize
513  END IF
514 
515  SELECT CASE(trim(key))
516  CASE('/timedisc/density')
517  sim%Timedisc%pvar%data4d(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX, &
518  sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%DENSITY) = 0.0
519 #ifndef PARALLEL
520  READ(unit) sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
521  sim%Mesh%KMIN:sim%Mesh%KMAX,sim%Physics%DENSITY)
522 #else
523  CALL mpi_file_set_view(handle, offset-1,default_mpi_real,filetype, 'native', mpi_info_null, ierror)
524  CALL mpi_file_read_all(handle, &
525  sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
526  sim%Mesh%KMIN:sim%Mesh%KMAX,sim%Physics%DENSITY), &
527  bufsize, default_mpi_real, status, ierror)
528  offset_0 = 0
529  CALL mpi_file_set_view(handle,offset_0,mpi_byte,mpi_byte,'native',mpi_info_null,ierror)
530 #endif
531  CASE('/timedisc/xvelocity')
532  sim%Timedisc%pvar%data4d(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX, &
533  sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%XVELOCITY) = 0.0
534 #ifndef PARALLEL
535  READ(unit) sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
536  sim%Mesh%KMIN:sim%Mesh%KMAX,sim%Physics%XVELOCITY)
537 #else
538  CALL mpi_file_set_view(handle, offset-1,default_mpi_real,filetype, 'native', mpi_info_null, ierror)
539  CALL mpi_file_read_all(handle, &
540  sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
541  sim%Mesh%KMIN:sim%Mesh%KMAX,sim%Physics%XVELOCITY),&
542  bufsize, default_mpi_real, status, ierror)
543  offset_0 = 0
544  CALL mpi_file_set_view(handle,offset_0,mpi_byte,mpi_byte,'native',mpi_info_null,ierror)
545 #endif
546  CASE('/timedisc/yvelocity')
547  sim%Timedisc%pvar%data4d(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX, &
548  sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%YVELOCITY) = 0.0
549 #ifndef PARALLEL
550  READ(unit) sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
551  sim%Mesh%KMIN:sim%Mesh%KMAX,sim%Physics%YVELOCITY)
552 #else
553  CALL mpi_file_set_view(handle, offset-1,default_mpi_real,filetype, 'native', mpi_info_null, ierror)
554  CALL mpi_file_read_all(handle, &
555  sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
556  sim%Mesh%KMIN:sim%Mesh%KMAX,sim%Physics%YVELOCITY),&
557  bufsize, default_mpi_real, status, ierror)
558  offset_0 = 0
559  CALL mpi_file_set_view(handle,offset_0,mpi_byte,mpi_byte,'native',mpi_info_null,ierror)
560 #endif
561  CASE('/timedisc/zvelocity')
562  sim%Timedisc%pvar%data4d(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX, &
563  sim%Mesh%KGMIN:sim%Mesh%KGMAX,sim%Physics%ZVELOCITY) = 0.0
564 #ifndef PARALLEL
565  READ(unit) sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
566  sim%Mesh%KMIN:sim%Mesh%KMAX,sim%Physics%ZVELOCITY)
567 #else
568  CALL mpi_file_set_view(handle, offset-1,default_mpi_real,filetype, 'native', mpi_info_null, ierror)
569  CALL mpi_file_read_all(handle, &
570  sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
571  sim%Mesh%KMIN:sim%Mesh%KMAX,sim%Physics%ZVELOCITY),&
572  bufsize, default_mpi_real, status, ierror)
573  offset_0 = 0
574  CALL mpi_file_set_view(handle,offset_0,mpi_byte,mpi_byte,'native',mpi_info_null,ierror)
575 #endif
576  CASE('/timedisc/pressure')
577  SELECT TYPE (phys => sim%Physics)
578  TYPE IS(physics_euler)
579  sim%Timedisc%pvar%data4d(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX, &
580  sim%Mesh%KGMIN:sim%Mesh%KGMAX,phys%PRESSURE) = 0.0
581 #ifndef PARALLEL
582  READ(unit) sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
583  sim%Mesh%KMIN:sim%Mesh%KMAX,phys%PRESSURE)
584 #else
585  CALL mpi_file_set_view(handle, offset-1,default_mpi_real,filetype, 'native', mpi_info_null, ierror)
586  CALL mpi_file_read_all(handle, &
587  sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
588  sim%Mesh%KMIN:sim%Mesh%KMAX,phys%PRESSURE),&
589  bufsize, default_mpi_real, status, ierror)
590  offset_0 = 0
591  CALL mpi_file_set_view(handle,offset_0,mpi_byte,mpi_byte,'native',mpi_info_null,ierror)
592 #endif
593  END SELECT
594  CASE('/physics/bccsound')
595  SELECT TYPE (phys => sim%Physics)
596  TYPE IS(physics_eulerisotherm)
597  phys%bccsound%data3d(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX, &
598  sim%Mesh%KGMIN:sim%Mesh%KGMAX) = 0.0
599 #ifndef PARALLEL
600  READ(unit) phys%bccsound%data3d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
601  sim%Mesh%KMIN:sim%Mesh%KMAX)
602 #else
603  CALL mpi_file_set_view(handle, offset-1,default_mpi_real,filetype, 'native', mpi_info_null, ierror)
604  CALL mpi_file_read_all(handle, &
605  phys%bccsound%data3d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
606  sim%Mesh%KMIN:sim%Mesh%KMAX),&
607  bufsize, default_mpi_real, status, ierror)
608  offset_0 = 0
609  CALL mpi_file_set_view(handle,offset_0,mpi_byte,mpi_byte,'native',mpi_info_null,ierror)
610 #endif
611  END SELECT
612  CASE('/physics/fcsound')
613  SELECT TYPE (phys => sim%Physics)
614  TYPE IS(physics_eulerisotherm)
615  phys%fcsound%data4d(sim%Mesh%IGMIN:sim%Mesh%IGMAX,sim%Mesh%JGMIN:sim%Mesh%JGMAX, &
616  sim%Mesh%KGMIN:sim%Mesh%KGMAX,1:sim%Mesh%NFACES) = 0.0
617 #ifndef PARALLEL
618  READ(unit) phys%fcsound%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
619  sim%Mesh%KMIN:sim%Mesh%KMAX,1:sim%Mesh%NFACES)
620 #else
621  CALL mpi_file_set_view(handle, offset-1,default_mpi_real,filetype, 'native', mpi_info_null, ierror)
622  CALL mpi_file_read_all(handle, &
623  phys%fcsound%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,sim%Mesh%JMIN:sim%Mesh%JMAX, &
624  sim%Mesh%KMIN:sim%Mesh%KMAX,1:sim%Mesh%NFACES),&
625  bufsize, default_mpi_real, status, ierror)
626  offset_0 = 0
627  CALL mpi_file_set_view(handle,offset_0,mpi_byte,mpi_byte,'native',mpi_info_null,ierror)
628 #endif
629  END SELECT
630  CASE('/sources/grav/binary/binpos')
631  srcptr => this%Sources
632  DO WHILE (ASSOCIATED(srcptr))
633  SELECT TYPE (gravity => srcptr)
634  TYPE IS (sources_gravity)
635  gravptr => gravity%glist
636  DO WHILE (ASSOCIATED(gravptr))
637  SELECT TYPE (binary => gravptr)
638  TYPE IS (gravity_binary)
639  READ(unit) binary%pos(1:3,1:2)
640  CLASS DEFAULT
641  ! do nothing or add fields/values in gravities
642  END SELECT
643  gravptr => gravptr%next
644  END DO
645  CLASS DEFAULT
646  ! do nothing or add fields/values in sources
647  END SELECT
648  srcptr => srcptr%next
649  END DO
650  CASE('/sources/grav/binary/mass')
651  srcptr => this%Sources
652  DO WHILE (ASSOCIATED(srcptr))
653  SELECT TYPE (gravity => srcptr)
654  TYPE IS (sources_gravity)
655  gravptr => gravity%glist
656  DO WHILE (ASSOCIATED(gravptr))
657  SELECT TYPE (binary => gravptr)
658  TYPE IS (gravity_binary)
659  READ(unit) binary%mass
660  CLASS DEFAULT
661  ! do nothing or add fields/values in gravities
662  END SELECT
663  gravptr => gravptr%next
664  END DO
665  CLASS DEFAULT
666  ! do nothing or add fields/values in sources
667  END SELECT
668  srcptr => srcptr%next
669  END DO
670  CASE('/sources/grav/binary/mass2')
671  srcptr => this%Sources
672  DO WHILE (ASSOCIATED(srcptr))
673  SELECT TYPE (gravity => srcptr)
674  TYPE IS (sources_gravity)
675  gravptr => gravity%glist
676  DO WHILE (ASSOCIATED(gravptr))
677  SELECT TYPE (binary => gravptr)
678  TYPE IS (gravity_binary)
679  READ(unit) binary%mass2
680  CLASS DEFAULT
681  ! do nothing or add fields/values in gravities
682  END SELECT
683  gravptr => gravptr%next
684  END DO
685  CLASS DEFAULT
686  ! do nothing or add fields/values in sources
687  END SELECT
688  srcptr => srcptr%next
689  END DO
690  CASE DEFAULT
691  SELECT CASE(l)
692  CASE(2)
693 #ifndef PARALLEL
694  ALLOCATE(ptr2(dims(1),dims(2)))
695  READ(unit) ptr2
696  DEALLOCATE(ptr2)
697 #endif
698  CASE(3)
699 #ifndef PARALLEL
700  ALLOCATE(ptr3(dims(1),dims(2),dims(3)))
701  READ(unit) ptr3
702  DEALLOCATE(ptr3)
703 #endif
704  CASE(4)
705 #ifndef PARALLEL
706  ALLOCATE(ptr4(dims(1),dims(2),dims(3),dims(4)))
707  READ(unit) ptr4
708  DEALLOCATE(ptr4)
709 #endif
710  CASE(5)
711 #ifndef PARALLEL
712  ALLOCATE(ptr5(dims(1),dims(2),dims(3),dims(4),dims(5)))
713  READ(unit) ptr5
714  DEALLOCATE(ptr5)
715 #endif
716  CASE DEFAULT
717  IF(bytes.GT.0) THEN
718 #ifndef PARALLEL
719  ALLOCATE(val(bytes))
720  READ(unit) val
721  DEALLOCATE(val)
722 #endif
723  END IF
724  END SELECT
725  END SELECT
726  offset = offset + bytes
727 
728 #ifndef PARALLEL
729  READ(unit, iostat=error) keylen
730 #else
731  buffer = ''
732  CALL mpi_file_seek(handle,offset-1,mpi_seek_set,ierror)
733  CALL mpi_file_read_all(handle,buffer(1:intsize),intsize,mpi_byte,status,ierror)
734  keylen = transfer(buffer(1:intsize),keylen)
735 
736  CALL mpi_file_get_position(handle,offset_0,ierror)
737  CALL mpi_file_get_size(handle,filesize,ierror)
738  IF (filesize.LE.offset_0) THEN
739  ierror=1
740  END IF
741  error = ierror
742 #endif
743  offset = offset + intsize
744  END DO
745 #ifndef PARALLEL
746  CLOSE(unit)
747 #else
748  CALL mpi_file_close(handle,ierror)
749 #endif
750 END SUBROUTINE loaddata
751 
752 END PROGRAM restart
integer, save default_mpi_real
default real type for MPI