PageRenderTime 7ms CodeModel.GetById 2ms app.highlight 110ms RepoModel.GetById 1ms app.codeStats 0ms

/input_output_hdf5.f90

https://bitbucket.org/kwmsmith/kaw-sim-hdf5
FORTRAN Modern | 819 lines | 647 code | 165 blank | 7 comment | 64 complexity | 5bef93028addc7bec2f9ba0f655e5e0a MD5 | raw file
  1module input_output_hdf5
  2    use globalpars
  3    use mpimod
  4    use io_constants, only : get_array_name
  5    use input_output_base_hdf5
  6    use params_io
  7    use table_funcs
  8    use HDF5
  9    use H5LT
 10    use H5TB
 11    implicit none
 12
 13    SAVE
 14
 15    PRIVATE
 16
 17    integer(HSIZE_T), dimension(1), parameter :: dummy_dims = (/ 20 /)
 18
 19    integer, parameter :: GRP_NAME_LEN = 20, NGPS = 17, NSTATS=11
 20
 21    integer, parameter :: ROOT = 0
 22
 23    character(LEN=GRP_NAME_LEN), dimension(NSTATS), parameter :: &
 24            stats_names = (/ "psi",  "vor",  "den", &
 25                             "phi",  "cur",  &
 26                             "bx",   "by",   &
 27                             "vx",   "vy",   &
 28                             "gdx",  "gdy" /)
 29
 30    character(LEN=GRP_NAME_LEN), dimension(NGPS), parameter :: &
 31            group_names = (/ "cpsi", "cvor", "cden", &
 32                     "psi",  "vor",  "den", &
 33                     "phi",  "cur",  &
 34                     "bx",   "by",   &
 35                     "vx",   "vy",   &
 36                     "gdx",  "gdy",  &
 37                     "scalars", "stats", &
 38                     "sim_params" /)
 39
 40    public :: init_hdf5, close_hdf5, create_file, init_file, write_rarr_hdf5,&
 41              write_carr_hdf5, read_carr_hdf5, open_file_hdf5, close_file_hdf5, HID_T, &
 42              dump_arrs_hdf5, load_cmplx_arrs_hdf5, group_names, DATA_FNAME, init_io_hdf5, &
 43              dump_scalars_hdf5, dump_params_hdf5, load_params_hdf5, input_params_t, &
 44              check_restart_params, find_last_arr_idx_hdf5
 45
 46    contains
 47
 48    subroutine check_restart_params(fname, input_params, error)
 49        character(*), intent(in) :: fname
 50        type(input_params_t), intent(inout) :: input_params
 51        integer, intent(out) :: error
 52        integer :: errbuf(1), ierr
 53
 54        type(input_params_t) :: hdf5_input_params
 55
 56        error = 0
 57
 58        if(pid .eq. 0) then
 59
 60            call load_params_hdf5(fname, hdf5_input_params, error)
 61
 62            error = 0
 63            ! we only allow NSTEPS to change, everything else must be the same.
 64            if(hdf5_input_params%nsteps .gt. input_params%nsteps ) then
 65                write(STDERR, *) "*** HDF5 NSTEPS > input params NSTEPS ***"
 66                call flush(STDERR)
 67                error = 1
 68            endif
 69            if(hdf5_input_params%max_wall_mins .ne. input_params%max_wall_mins ) then
 70                write(STDERR, *) "*** HDF5 param does not equal input param ***"
 71                call flush(STDERR)
 72                error = 1
 73                            endif
 74            if(hdf5_input_params%nout_arrs .ne. input_params%nout_arrs ) then
 75                write(STDERR, *) "*** HDF5 param does not equal input param ***"
 76                call flush(STDERR)
 77                error = 1
 78                            endif
 79            if(hdf5_input_params%nout_scals .ne. input_params%nout_scals ) then
 80                write(STDERR, *) "*** HDF5 param does not equal input param ***"
 81                call flush(STDERR)
 82                error = 1
 83                            endif
 84            if(hdf5_input_params%rng_seed .ne. input_params%rng_seed ) then
 85                write(STDERR, *) "*** HDF5 param does not equal input param ***"
 86                call flush(STDERR)
 87                error = 1
 88                            endif
 89
 90            if(hdf5_input_params%tstep .ne. input_params%tstep ) then
 91                write(STDERR, *) "*** HDF5 param does not equal input param ***"
 92                call flush(STDERR)
 93                error = 1
 94                            endif
 95            if(hdf5_input_params%lx .ne. input_params%lx ) then
 96                write(STDERR, *) "*** HDF5 param does not equal input param ***"
 97                call flush(STDERR)
 98                error = 1
 99                            endif
100            if(hdf5_input_params%ly .ne. input_params%ly ) then
101                write(STDERR, *) "*** HDF5 param does not equal input param ***"
102                call flush(STDERR)
103                error = 1
104                            endif
105            if(hdf5_input_params%viscos .ne. input_params%viscos ) then
106                write(STDERR, *) "*** HDF5 param does not equal input param ***"
107                call flush(STDERR)
108                error = 1
109                            endif
110            if(hdf5_input_params%h_viscos .ne. input_params%h_viscos ) then
111                write(STDERR, *) "*** HDF5 param does not equal input param ***"
112                call flush(STDERR)
113                error = 1
114                            endif
115            if(hdf5_input_params%resis .ne. input_params%resis ) then
116                write(STDERR, *) "*** HDF5 param does not equal input param ***"
117                call flush(STDERR)
118                error = 1
119                            endif
120            if(hdf5_input_params%h_resis .ne. input_params%h_resis ) then
121                write(STDERR, *) "*** HDF5 param does not equal input param ***"
122                call flush(STDERR)
123                error = 1
124                            endif
125            if(hdf5_input_params%diffus .ne. input_params%diffus ) then
126                write(STDERR, *) "*** HDF5 param does not equal input param ***"
127                call flush(STDERR)
128                error = 1
129                            endif
130            if(hdf5_input_params%h_diffus .ne. input_params%h_diffus ) then
131                write(STDERR, *) "*** HDF5 param does not equal input param ***"
132                call flush(STDERR)
133                error = 1
134                            endif
135            if(hdf5_input_params%spectrum_slope .ne. input_params%spectrum_slope) then
136                write(STDERR, *) "*** HDF5 param does not equal input param ***"
137                call flush(STDERR)
138                error = 1
139                            endif
140            if(hdf5_input_params%spectrum_peak .ne. input_params%spectrum_peak ) then
141                write(STDERR, *) "*** HDF5 param does not equal input param ***"
142                call flush(STDERR)
143                error = 1
144                            endif
145            if(hdf5_input_params%eb .ne. input_params%eb ) then
146                write(STDERR, *) "*** HDF5 param does not equal input param ***"
147                call flush(STDERR)
148                error = 1
149                            endif
150            if(hdf5_input_params%ev .ne. input_params%ev ) then
151                write(STDERR, *) "*** HDF5 param does not equal input param ***"
152                call flush(STDERR)
153                error = 1
154                            endif
155            if(hdf5_input_params%en .ne. input_params%en ) then
156                write(STDERR, *) "*** HDF5 param does not equal input param ***"
157                call flush(STDERR)
158                error = 1
159                            endif
160            if(hdf5_input_params%rho_s2 .ne. input_params%rho_s2 ) then
161                write(STDERR, *) "*** HDF5 param does not equal input param ***"
162                call flush(STDERR)
163                error = 1
164                            endif
165            if(hdf5_input_params%dtfac .ne. input_params%dtfac ) then
166                write(STDERR, *) "*** HDF5 param does not equal input param ***"
167                call flush(STDERR)
168                error = 1
169                            endif
170            if(hdf5_input_params%mason .ne. input_params%mason ) then
171                write(STDERR, *) "*** HDF5 param does not equal input param ***"
172                call flush(STDERR)
173                error = 1
174                            endif
175            if(hdf5_input_params%famp .ne. input_params%famp ) then
176                write(STDERR, *) "*** HDF5 param does not equal input param ***"
177                call flush(STDERR)
178                error = 1
179                            endif
180                errbuf(1) = error
181        endif
182
183        CALL MPI_BCAST(errbuf, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
184        error = errbuf(1)
185
186    end subroutine check_restart_params
187
188
189    subroutine find_last_arr_idx_hdf5(fname, last_idx, error)
190        implicit none
191        character(*), intent(in) :: fname
192        integer, intent(out) :: last_idx
193        integer, intent(out) :: error
194
195        integer(HID_T) :: file_id, group_id
196
197        integer :: exists, nout_arrs, i, ierr, intbuf(1)
198
199        if(pid .eq. 0) then
200
201            call open_file_hdf5(fname, file_id, error)
202
203            call h5gopen_f(file_id, "sim_params", group_id, error)
204            call read_attr_integer(group_id, "NOUT_ARRS", nout_arrs, error)
205            call h5gclose_f(group_id, error)
206
207            call h5gopen_f(file_id, "cden", group_id, error)
208
209            last_idx = -1
210            i = 0
211            do
212                exists = h5ltfind_dataset_f(group_id, get_array_name(i))
213                if (exists .eq. 1) then
214                    last_idx = i
215                else
216                    exit
217                endif
218                i = i + nout_arrs
219            enddo
220
221            call h5gclose_f(group_id, error)
222
223            call close_file_hdf5(file_id, error)
224
225            intbuf(1) = last_idx
226
227        endif
228
229        CALL MPI_BCAST(intbuf, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
230        last_idx = intbuf(1)
231
232    end subroutine find_last_arr_idx_hdf5
233
234    subroutine load_params_hdf5(fname, input_params, error)
235        implicit none
236        character(*), intent(in) :: fname
237        type(input_params_t), intent(out) :: input_params
238        integer, intent(out) :: error
239
240        integer(HID_T) :: file_id, group_id, input_int
241
242        if(pid .ne. 0) return
243
244        call open_file_hdf5(fname, file_id, error)
245        call h5gopen_f(file_id, "sim_params", group_id, error)
246
247        call read_attr_integer(group_id, "MAJOR_VERSION", input_int, error)
248        if (input_int .ne. MAJOR_VERSION) then
249            write(STDERR, *) "*** datafile MAJOR_VERSION ", &
250                             input_int, " .ne. to compiled MAJOR_VERSION ", &
251                             MAJOR_VERSION, " ***"
252            call flush(STDERR)
253            error = 1
254        endif
255
256        call read_attr_integer(group_id, "MINOR_VERSION", input_int, error)
257        if (input_int .ne. MINOR_VERSION) then
258            write(STDERR, *) "*** datafile MINOR_VERSION ", &
259                             input_int, " .ne. to compiled MINOR_VERSION ", &
260                             MINOR_VERSION, " ***"
261            call flush(STDERR)
262            error = 1
263        endif
264
265        call read_attr_integer(group_id, "NX", input_int, error)
266        if (input_int .ne. NX) then
267            write(STDERR, *) "*** datafile NX ", &
268                             input_int, " .ne. to compiled NX ", &
269                             NX, " ***"
270            call flush(STDERR)
271            error = 1
272        endif
273
274        call read_attr_integer(group_id, "NY", input_int, error)
275        if (input_int .ne. NY) then
276            write(STDERR, *) "*** datafile NY ", &
277                             input_int, " .ne. to compiled NY ", &
278                             NY, " ***"
279            call flush(STDERR)
280            error = 1
281        endif
282
283        call read_attr_integer(group_id, "NP", input_int, error)
284        if (input_int .ne. NP) then
285            write(STDERR, *) "*** datafile NP ", &
286                             input_int, " .ne. to compiled NP ", &
287                             NP, " ***"
288            call flush(STDERR)
289            error = 1
290        endif
291
292        call read_attr_integer(group_id, "NSTEPS", input_params%nsteps, error)
293        call read_attr_integer(group_id, "MAX_WALL_MINS", input_params%max_wall_mins, error)
294        call read_attr_integer(group_id, "NOUT_ARRS", input_params%nout_arrs, error)
295        call read_attr_integer(group_id, "NOUT_SCALS", input_params%nout_scals, error)
296        call read_attr_integer(group_id, "RNG_SEED", input_params%rng_seed, error)
297
298        call read_attr_real(group_id, "TSTEP", input_params%tstep, error)
299        call read_attr_real(group_id, "LX", input_params%lx, error)
300        call read_attr_real(group_id, "LY", input_params%ly, error)
301        call read_attr_real(group_id, "VISCOS", input_params%viscos, error)
302        call read_attr_real(group_id, "H_VISCOS", input_params%h_viscos, error)
303        call read_attr_real(group_id, "RESIS", input_params%resis, error)
304        call read_attr_real(group_id, "H_RESIS", input_params%h_resis, error)
305        call read_attr_real(group_id, "DIFFUS", input_params%diffus, error)
306        call read_attr_real(group_id, "H_DIFFUS", input_params%h_diffus, error)
307        call read_attr_real(group_id, "SPECTRUM_SLOPE", input_params%spectrum_slope, error)
308        call read_attr_real(group_id, "SPECTRUM_PEAK", input_params%spectrum_peak, error)
309        call read_attr_real(group_id, "EB", input_params%eb, error)
310        call read_attr_real(group_id, "EV", input_params%ev, error)
311        call read_attr_real(group_id, "EN", input_params%en, error)
312        call read_attr_real(group_id, "RHO_S2", input_params%rho_s2, error)
313        call read_attr_real(group_id, "DTFAC", input_params%dtfac, error)
314        call read_attr_real(group_id, "MASON", input_params%mason, error)
315        call read_attr_real(group_id, "FAMP", input_params%famp, error)
316
317        call h5gclose_f(group_id, error)
318        call close_file_hdf5(file_id, error)
319
320    end subroutine load_params_hdf5
321
322    subroutine read_attr_real(group_id, attr_name, attr_value, err)
323        integer(HID_T), intent(in) :: group_id
324        character(*), intent(in) :: attr_name
325        real(sp), intent(out) :: attr_value
326        integer, intent(out) :: err
327
328        integer(HID_T) :: attr_id
329
330        call h5aopen_f(group_id, attr_name, attr_id, err)
331        call h5aread_f(attr_id, H5T_NATIVE_REAL, attr_value, &
332                        dummy_dims, err)
333        call h5aclose_f(attr_id, err)
334    end subroutine read_attr_real
335
336    subroutine read_attr_integer(group_id, attr_name, attr_value, err)
337        integer(HID_T), intent(in) :: group_id
338        character(*), intent(in) :: attr_name
339        integer, intent(out) :: attr_value
340        integer, intent(out) :: err
341
342        integer(HID_T) :: attr_id
343
344        call h5aopen_f(group_id, attr_name, attr_id, err)
345        call h5aread_f(attr_id, H5T_NATIVE_INTEGER, attr_value, &
346                        dummy_dims, err)
347        call h5aclose_f(attr_id, err)
348
349    end subroutine read_attr_integer
350
351    subroutine read_attr_character(group_id, attr_name, attr_value, err)
352        integer(HID_T), intent(in) :: group_id
353        character(*), intent(in) :: attr_name
354        character(*), intent(out) :: attr_value
355        integer, intent(out) :: err
356
357        character(STR_LEN) :: buf
358
359        integer(HID_T) :: attr_id, atype_id
360
361        call h5tcopy_f(H5T_NATIVE_CHARACTER, atype_id, err)
362        call h5tset_size_f(atype_id, int(STR_LEN, kind=HSIZE_T), err)
363
364        call h5aopen_f(group_id, attr_name, attr_id, err)
365        call h5aread_f(attr_id, atype_id, buf, &
366                        dummy_dims, err)
367        call h5aclose_f(attr_id, err)
368
369        attr_value = buf
370
371    end subroutine read_attr_character
372
373    subroutine dump_params_restart(fname, input_params)
374        character(*), intent(in) :: fname
375        type(input_params_t), intent(in) :: input_params
376
377        integer(HID_T) :: file_id, group_id
378
379        integer :: err
380
381        if(pid .ne. 0) return
382
383        call open_file_hdf5(fname, file_id, err)
384        call h5gopen_f(file_id, "sim_params", group_id, err)
385
386        call write_attr_integer(group_id, "NSTEPS", input_params%nsteps, err)
387
388        call h5gclose_f(group_id, err)
389        call close_file_hdf5(file_id, err)
390
391        contains
392
393        subroutine write_attr_integer(group_id, attr_name, attr_value, err)
394            integer(HID_T), intent(in) :: group_id
395            character(*), intent(in) :: attr_name
396            integer, intent(in) :: attr_value
397            integer, intent(out) :: err
398
399            integer(HID_T) :: attr_id
400
401            call h5aopen_f(group_id, attr_name, attr_id, err)
402            call h5awrite_f(attr_id, H5T_NATIVE_INTEGER,&
403                            attr_value, dummy_dims, err)
404            call h5aclose_f(attr_id, err)
405
406        end subroutine write_attr_integer
407
408    end subroutine dump_params_restart
409
410    subroutine dump_params_hdf5(fname, input_params, restart)
411        implicit none
412        character(*), intent(in) :: fname
413        type(input_params_t), intent(in) :: input_params
414        logical, intent(in) :: restart
415
416        integer(HID_T) :: atype_id, space_id, file_id, group_id
417        integer :: err
418
419        integer(HSIZE_T), dimension(1), parameter :: dummy_dims = (/ 20 /)
420
421        if(pid .ne. 0) return
422
423        if(restart) then
424            call dump_params_restart(fname, input_params)
425            return
426        endif
427
428        call open_file_hdf5(fname, file_id, err)
429        call h5gopen_f(file_id, "sim_params", group_id, err)
430
431        call h5screate_f(H5S_SCALAR_F, space_id, err)
432        CALL h5tcopy_f(H5T_NATIVE_CHARACTER, atype_id, err)
433        CALL h5tset_size_f(atype_id, int(STR_LEN, kind=HSIZE_T), err)
434
435        call write_attr_integer("MAJOR_VERSION", MAJOR_VERSION, err)
436        call write_attr_integer("MINOR_VERSION", MINOR_VERSION, err)
437        call write_attr_integer("NX", nx, err)
438        call write_attr_integer("NY", ny, err)
439        call write_attr_integer("NP", np, err)
440
441        ! REVISION, INTEGRATOR & MODEL -- character arrays
442        call write_attr_character("REVISION", REVISION, err)
443        call write_attr_character("INTEGRATOR", INTEGRATOR_, err)
444        call write_attr_character("MODEL", MODEL, err)
445
446        call write_attr_integer("NSTEPS", input_params%nsteps, err)
447        call write_attr_integer("MAX_WALL_MINS", input_params%max_wall_mins, err)
448        call write_attr_integer("NOUT_ARRS", input_params%nout_arrs, err)
449        call write_attr_integer("NOUT_SCALS", input_params%nout_scals, err)
450        call write_attr_integer("RNG_SEED", input_params%rng_seed, err)
451
452        call write_attr_real("TSTEP", input_params%tstep, err)
453        call write_attr_real("LX", input_params%lx, err)
454        call write_attr_real("LY", input_params%ly, err)
455        call write_attr_real("VISCOS", input_params%viscos, err)
456        call write_attr_real("H_VISCOS", input_params%h_viscos, err)
457        call write_attr_real("RESIS", input_params%resis, err)
458        call write_attr_real("H_RESIS", input_params%h_resis, err)
459        call write_attr_real("DIFFUS", input_params%diffus, err)
460        call write_attr_real("H_DIFFUS", input_params%h_diffus, err)
461        call write_attr_real("SPECTRUM_SLOPE", input_params%spectrum_slope, err)
462        call write_attr_real("SPECTRUM_PEAK", input_params%spectrum_peak, err)
463        call write_attr_real("EB", input_params%eb, err)
464        call write_attr_real("EV", input_params%ev, err)
465        call write_attr_real("EN", input_params%en, err)
466        call write_attr_real("RHO_S2", input_params%rho_s2, err)
467        call write_attr_real("DTFAC", input_params%dtfac, err)
468        call write_attr_real("MASON", input_params%mason, err)
469        call write_attr_real("FAMP", input_params%famp, err)
470
471        call h5sclose_f(space_id, err)
472        call h5gclose_f(group_id, err)
473        call close_file_hdf5(file_id, err)
474
475        contains!{{{
476
477        subroutine write_attr_real(attr_name, attr_value, err)
478            character(*), intent(in) :: attr_name
479            real(sp), intent(in) :: attr_value
480            integer, intent(out) :: err
481
482            integer(HID_T) :: attr_id
483
484            call h5acreate_f(group_id, attr_name,&
485                             H5T_NATIVE_REAL, space_id, attr_id, err)
486            call h5awrite_f(attr_id, H5T_NATIVE_REAL,&
487                            attr_value, dummy_dims, err)
488            call h5aclose_f(attr_id, err)
489        end subroutine write_attr_real
490
491        subroutine write_attr_integer(attr_name, attr_value, err)
492            character(*), intent(in) :: attr_name
493            integer, intent(in) :: attr_value
494            integer, intent(out) :: err
495
496            integer(HID_T) :: attr_id
497
498            call h5acreate_f(group_id, attr_name,&
499                             H5T_NATIVE_INTEGER, space_id, attr_id, err)
500            call h5awrite_f(attr_id, H5T_NATIVE_INTEGER,&
501                            attr_value, dummy_dims, err)
502            call h5aclose_f(attr_id, err)
503        end subroutine write_attr_integer
504
505        subroutine write_attr_character(attr_name, attr_value, err)
506            character(*), intent(in) :: attr_name
507            character(*), intent(in) :: attr_value
508            integer, intent(out) :: err
509
510            character(STR_LEN) :: buf
511
512            integer(HID_T) :: attr_id
513
514            buf = attr_value
515        
516            call h5acreate_f(group_id, attr_name,&
517                             atype_id, space_id, attr_id, err)
518            call h5awrite_f(attr_id, atype_id,&
519                            buf, dummy_dims, err)
520            call h5aclose_f(attr_id, err)
521        end subroutine write_attr_character
522!}}}
523
524    end subroutine dump_params_hdf5
525
526    subroutine init_io_hdf5(fname, overwrite, restart, error)
527        implicit none
528        character(*), intent(in) :: fname
529        logical, intent(in) :: restart
530        logical, intent(in) :: overwrite
531        integer, intent(out) :: error
532        integer(HID_T) :: file_id, group_id
533        integer :: i
534
535        call init_table_funcs()
536
537        if(restart) return
538
539        call create_file(fname, overwrite, error)
540        if(error .ne. 0) return
541        call init_file(fname, group_names, error)
542
543        if(pid .eq. ROOT) then
544            call open_file_hdf5(fname, file_id, error)
545            call h5gopen_f(file_id, "scalars", group_id, error)
546            call make_table_scalars(group_id, error)
547            call h5gclose_f(group_id, error)
548            call h5gopen_f(file_id, "stats", group_id, error)
549            do i = 1, NSTATS
550                call make_table_stats(group_id, stats_names(i), error)
551            enddo
552            call h5gclose_f(group_id, error)
553            call close_file_hdf5(file_id, error)
554        endif
555
556    end subroutine init_io_hdf5
557
558    subroutine close_io_hdf5()
559        implicit none
560    end subroutine close_io_hdf5
561
562    subroutine dump_arrs_hdf5(fname, step_num, nout_arrs, cfields, iserr)
563        implicit none
564        character(*), intent(in) :: fname
565        integer, intent(in) :: step_num, nout_arrs
566        complex(sp), dimension(cdim, NF), intent(in) :: cfields
567        integer, intent(out) :: iserr
568
569        iserr = 0
570        if(nout_arrs .eq. 0) then
571            write(STDERR, *) "*** Error: nout_arrs is 0. ***"
572            call flush(STDERR)
573            iserr = 1
574            return
575        endif
576
577        if(mod(step_num, nout_arrs) .eq. 0) then
578            call dump_cmplx_arrs(fname, step_num, cfields)
579            call dump_real_arrs(fname, step_num, cfields)
580        endif
581
582    end subroutine dump_arrs_hdf5
583
584    subroutine load_cmplx_arrs_hdf5(fname, snap, cfield, error)
585        implicit none
586        character(*), intent(in) :: fname
587        integer, intent(in) :: snap
588        complex(sp), dimension(cdim, NF), intent(out) :: cfield
589        integer, intent(out) :: error
590
591        integer(HID_T) :: file_id
592
593        error = 0
594        call open_file_hdf5(fname, file_id, error)
595        if (error .ne. 0) return
596
597        call read_carr_hdf5(file_id, "cpsi", snap, cfield(:,PSI_IDX), error)
598        if(error .ne. 0) goto 90
599
600        call read_carr_hdf5(file_id, "cvor", snap, cfield(:,VOR_IDX), error)
601        if(error .ne. 0) goto 90
602
603        call read_carr_hdf5(file_id, "cden", snap, cfield(:,DEN_IDX), error)
604        if(error .ne. 0) goto 90
605
606        90 call close_file_hdf5(file_id, error)
607
608    end subroutine load_cmplx_arrs_hdf5
609
610    subroutine dump_cmplx_arrs(fname, snap, cfield)
611        implicit none
612        character(*), intent(in) :: fname
613        integer, intent(in) :: snap
614        complex(sp), dimension(cdim, NF), intent(in) :: cfield
615
616        integer(HID_T) :: file_id
617        integer :: error
618
619        error = 0
620        call open_file_hdf5(fname, file_id, error)
621        if(error .ne. 0) return
622
623        call write_carr_hdf5(file_id, "cpsi", snap, cfield(:,PSI_IDX), error)
624        if(error .ne. 0) goto 70
625
626        call write_carr_hdf5(file_id, "cvor", snap, cfield(:,VOR_IDX), error)
627        if(error .ne. 0) goto 70
628
629        call write_carr_hdf5(file_id, "cden", snap, cfield(:,DEN_IDX), error)
630        if(error .ne. 0) goto 70
631
632        70 call close_file_hdf5(file_id, error)
633    end subroutine dump_cmplx_arrs
634
635    subroutine dump_real_arrs(fname, snap, cfield)
636        use pfft
637        use spectral
638        implicit none
639        character(*), intent(in) :: fname
640        integer, intent(in) :: snap
641        complex(sp), dimension(cdim, NF), intent(in) :: cfield
642
643        integer(HID_T) :: file_id
644        integer :: error
645        real(sp), dimension(ldim) :: r_aux
646
647        error = 0
648        call open_file_hdf5(fname, file_id, error)
649        if(error .ne. 0) return
650
651        call pfft_c2r(cfield(:,PSI_IDX), r_aux)
652        call write_rarr_hdf5(file_id, "psi", snap, r_aux, error)
653        if (error .ne. 0) goto 80
654
655        call pfft_c2r(cfield(:,VOR_IDX), r_aux)
656        call write_rarr_hdf5(file_id, "vor", snap, r_aux, error)
657        if (error .ne. 0) goto 80
658
659        call pfft_c2r(cfield(:,DEN_IDX), r_aux)
660        call write_rarr_hdf5(file_id, "den", snap, r_aux, error)
661        if (error .ne. 0) goto 80
662
663        call pfft_c2r(-km2*cfield(:,VOR_IDX), r_aux)
664        call write_rarr_hdf5(file_id, "phi", snap, r_aux, error)
665        if (error .ne. 0) goto 80
666
667        call pfft_c2r(-k2*cfield(:,PSI_IDX), r_aux)
668        call write_rarr_hdf5(file_id, "cur", snap, r_aux, error)
669        if (error .ne. 0) goto 80
670
671        call derivative(X_DIR, cfield(:,PSI_IDX), r_aux)
672        call write_rarr_hdf5(file_id, "bx", snap, r_aux, error)
673        if (error .ne. 0) goto 80
674
675        call derivative(Y_DIR, cfield(:,PSI_IDX), r_aux)
676        call write_rarr_hdf5(file_id, "by", snap, r_aux, error)
677        if (error .ne. 0) goto 80
678
679        call derivative(X_DIR, cfield(:,DEN_IDX), r_aux)
680        call write_rarr_hdf5(file_id, "gdx", snap, r_aux, error)
681        if (error .ne. 0) goto 80
682
683        call derivative(Y_DIR, cfield(:,DEN_IDX), r_aux)
684        call write_rarr_hdf5(file_id, "gdy", snap, r_aux, error)
685        if (error .ne. 0) goto 80
686
687        call derivative(X_DIR, -km2*cfield(:,VOR_IDX), r_aux)
688        call write_rarr_hdf5(file_id, "vx", snap, r_aux, error)
689        if (error .ne. 0) goto 80
690
691        call derivative(Y_DIR, -km2*cfield(:,VOR_IDX), r_aux)
692        call write_rarr_hdf5(file_id, "vy", snap, r_aux, error)
693        if (error .ne. 0) goto 80
694
695        80 call close_file_hdf5(file_id, error)
696    end subroutine dump_real_arrs
697
698    subroutine dump_stats_hdf5(fname, dump_idx, cfields)
699        use spectral
700        use pfft
701        use statistics
702        implicit none
703        character(*), intent(in) :: fname
704        integer, intent(in) :: dump_idx
705        complex(sp), dimension(cdim, NF), intent(in) :: cfields
706
707        integer :: error
708        integer(HID_T) :: file_id, group_id
709
710        if(pid .eq. ROOT) then
711            error = 0
712            call open_file_hdf5(fname, file_id, error)
713            call h5gopen_f(file_id, "stats", group_id, error)
714        endif
715
716        call write_single_("psi", cfields(:,PSI_IDX))
717        call write_single_("phi", -km2*cfields(:,VOR_IDX))
718        call write_single_("cur", -k2*cfields(:,PSI_IDX))
719
720        call write_single_("bx", ikx*cfields(:,PSI_IDX))
721        call write_single_("by", iky*cfields(:,PSI_IDX))
722
723        call write_single_("den", cfields(:,DEN_IDX))
724        call write_single_("gdx", ikx*cfields(:,DEN_IDX))
725        call write_single_("gdy", iky*cfields(:,DEN_IDX))
726
727        call write_single_("vor", cfields(:,VOR_IDX))
728        call write_single_("vx", -km2*ikx*cfields(:,VOR_IDX))
729        call write_single_("vy", -km2*iky*cfields(:,VOR_IDX))
730
731        if(pid .eq. ROOT) then
732            call h5gclose_f(group_id, error)
733            call close_file_hdf5(file_id, error)
734        endif
735
736        contains
737
738        subroutine write_single_(nm, carr)
739            implicit none
740            character(*), intent(in) :: nm
741            complex(sp), dimension(cdim), intent(in) :: carr
742            real(sp), dimension(ldim) :: rfield
743            real(db) :: ave, std_dev, skew, kurt
744            integer :: nrecords
745            type(stat_record), dimension(1) :: data_
746
747            call pfft_c2r(carr, rfield)
748            call get_moments(rfield, ave, std_dev, skew, kurt)
749
750            if(pid .eq. ROOT) then
751                error = 0
752                nrecords = 1
753                data_ = (/ stat_record(dump_idx, ave, std_dev, skew, kurt) /)
754                call append_record_stats(group_id, nm, nrecords, data_, error)
755            endif
756
757        end subroutine write_single_
758
759    end subroutine dump_stats_hdf5
760
761    subroutine dump_scalars_hdf5(fname, dump_idx, time, basic_steps, rho_s2,&
762                                 cfield)
763        implicit none
764        character(*), intent(in) :: fname
765        integer, intent(in) :: dump_idx, basic_steps
766        real(db), intent(in) :: time
767        real(sp), intent(in) :: rho_s2
768        complex(sp), dimension(cdim, NF), intent(in) :: cfield
769
770        call dump_scalars_hdf5_(fname, dump_idx, time, basic_steps, rho_s2,&
771                                cfield)
772        call dump_stats_hdf5(fname, dump_idx, cfield)
773    end subroutine dump_scalars_hdf5
774
775    subroutine dump_scalars_hdf5_(fname, dump_idx, time, basic_steps, rho_s2,&
776                                cfield)
777        use spectral
778        implicit none
779        character(*), intent(in) :: fname
780        integer, intent(in) :: dump_idx, basic_steps
781        real(db), intent(in) :: time
782        real(sp), intent(in) :: rho_s2
783        complex(sp), dimension(cdim, NF), intent(in) :: cfield
784        real(db) :: msf, be, ve, ne
785        integer :: root, error, nrecords
786        type(scalar_record), dimension(1) :: data_
787
788        integer(HID_T) :: file_id, group_id
789
790        root = 0
791
792        ! msf
793        call reduce_sum(root, cfield(:,PSI_IDX)*conjg(cfield(:,PSI_IDX)), msf)
794
795        ! b-energy
796        call reduce_sum(root, k2*cfield(:,PSI_IDX)*conjg(cfield(:,PSI_IDX)), be)
797
798        ! v-energy
799        call reduce_sum(root, km2*cfield(:,VOR_IDX)*conjg(cfield(:,VOR_IDX)), ve)
800
801        ! internal-energy
802        call reduce_sum(root, cfield(:,DEN_IDX)*conjg(cfield(:,DEN_IDX)), ne)
803        ne = rho_s2 * ne
804
805
806        if(pid .eq. root) then
807            error = 0
808            nrecords = 1
809            data_ = (/ scalar_record(dump_idx, mpi_wtime()-start_time,&
810                                     basic_steps, time, be, ve, ne, msf) /)
811            call open_file_hdf5(fname, file_id, error)
812            call h5gopen_f(file_id, "scalars", group_id, error)
813            call append_records_scalars(group_id, nrecords, data_, error)
814            call h5gclose_f(group_id, error)
815            call close_file_hdf5(file_id, error)
816        endif
817    end subroutine dump_scalars_hdf5_
818
819end module input_output_hdf5