All Classes Files Functions Variables Pages
module_snapshot_surface.f90
Go to the documentation of this file.
1 !=====================================================================================================================================
2 ! EFISPEC3D !
3 ! (Elements FInis SPECtraux 3D) !
4 ! !
5 ! http://efispec.free.fr !
6 ! !
7 ! !
8 ! This file is part of EFISPEC3D !
9 ! Please refer to http://efispec.free.fr if you use it or part of it !
10 ! !
11 ! !
12 !1 ---> French License: CeCILL V2 !
13 ! !
14 ! Copyright BRGM 2009 contributeurs : Florent DE MARTIN !
15 ! David MICHEA !
16 ! Philippe THIERRY !
17 ! !
18 ! Contact: f.demartin at brgm.fr !
19 ! !
20 ! Ce logiciel est un programme informatique servant a resoudre l'equation du !
21 ! mouvement en trois dimensions via une methode des elements finis spectraux. !
22 ! !
23 ! Ce logiciel est regi par la licence CeCILL soumise au droit francais et !
24 ! respectant les principes de diffusion des logiciels libres. Vous pouvez !
25 ! utiliser, modifier et/ou redistribuer ce programme sous les conditions de la !
26 ! licence CeCILL telle que diffusee par le CEA, le CNRS et l'INRIA sur le site !
27 ! "http://www.cecill.info". !
28 ! !
29 ! En contrepartie de l'accessibilite au code source et des droits de copie, de !
30 ! modification et de redistribution accordes par cette licence, il n'est offert !
31 ! aux utilisateurs qu'une garantie limitee. Pour les memes raisons, seule une !
32 ! responsabilite restreinte pese sur l'auteur du programme, le titulaire des !
33 ! droits patrimoniaux et les concedants successifs. !
34 ! !
35 ! A cet egard l'attention de l'utilisateur est attiree sur les risques associes !
36 ! au chargement, a l'utilisation, a la modification et/ou au developpement et a !
37 ! la reproduction du logiciel par l'utilisateur etant donne sa specificite de !
38 ! logiciel libre, qui peut le rendre complexe a manipuler et qui le reserve donc !
39 ! a des developpeurs et des professionnels avertis possedant des connaissances !
40 ! informatiques approfondies. Les utilisateurs sont donc invites a charger et !
41 ! tester l'adequation du logiciel a leurs besoins dans des conditions permettant !
42 ! d'assurer la securite de leurs systemes et ou de leurs donnees et, plus !
43 ! generalement, a l'utiliser et l'exploiter dans les memes conditions de !
44 ! securite. !
45 ! !
46 ! Le fait que vous puissiez acceder a cet en-tete signifie que vous avez pris !
47 ! connaissance de la licence CeCILL et que vous en avez accepte les termes. !
48 ! !
49 ! !
50 !2 ---> International license: GNU GPL V3 !
51 ! !
52 ! EFISPEC3D is a computer program that solves the three-dimensional equations of !
53 ! motion using a finite spectral-element method. !
54 ! !
55 ! Copyright (C) 2009 Florent DE MARTIN !
56 ! !
57 ! Contact: f.demartin at brgm.fr !
58 ! !
59 ! This program is free software: you can redistribute it and/or modify it under !
60 ! the terms of the GNU General Public License as published by the Free Software !
61 ! Foundation, either version 3 of the License, or (at your option) any later !
62 ! version. !
63 ! !
64 ! This program is distributed in the hope that it will be useful, but WITHOUT ANY !
65 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A !
66 ! PARTICULAR PURPOSE. See the GNU General Public License for more details. !
67 ! !
68 ! You should have received a copy of the GNU General Public License along with !
69 ! this program. If not, see http://www.gnu.org/licenses/. !
70 ! !
71 ! !
72 !3 ---> Third party libraries !
73 ! !
74 ! EFISPEC3D uses the following of third party libraries: !
75 ! !
76 ! --> METIS 5.1.0 !
77 ! see http://glaros.dtc.umn.edu/gkhome/metis/metis/overview !
78 ! !
79 ! --> Lib_VTK_IO !
80 ! see S. Zaghi's website: https://github.com/szaghi/Lib_VTK_IO !
81 ! !
82 ! --> INTERP_LINEAR !
83 ! see J. Burkardt website: http://people.sc.fsu.edu/~jburkardt/ !
84 ! !
85 ! --> SAC !
86 ! http://ds.iris.edu/files/sac-manual/ !
87 ! !
88 ! --> EXODUS II !
89 ! http://sourceforge.net/projects/exodusii/ !
90 ! !
91 ! --> NETCDF !
92 ! http://www.unidata.ucar.edu/software/netcdf/ !
93 ! !
94 ! --> HDF5 !
95 ! http://www.hdfgroup.org/HDF5/ !
96 ! !
97 ! Some of these libraries are located in directory lib and pre-compiled !
98 ! with intel compiler for x86_64 platform. !
99 ! !
100 ! !
101 !4 ---> Related Articles !
102 ! !
103 ! De Martin, F., Matsushima, M., Kawase, H. (BSSA, 2013) !
104 ! Impact of geometric effects on near-surface Green's functions !
105 ! doi:10.1785/0120130039 !
106 ! !
107 ! Aochi, H., Ducellier, A., Dupros, F., Delatre, M., Ulrich, T., De Martin, F., !
108 ! Yoshimi, M., (Pure Appl. Geophys. 2013) !
109 ! Finite Difference Simulations of Seismic Wave Propagation for the 2007 Mw 6.6 !
110 ! Niigata-ken Chuetsu-Oki Earthquake: Validity of Models and Reliable !
111 ! Input Ground Motion in the Near-Field !
112 ! doi:10.1007/s00024-011-0429-5 !
113 ! !
114 ! De Martin, F. (BSSA, 2011) !
115 ! Verification of a Spectral-Element Method Code for the Southern California !
116 ! Earthquake Center LOH.3 Viscoelastic Case !
117 ! doi:10.1785/0120100305 !
118 ! !
119 !=====================================================================================================================================
120 
123 
126 module mod_snapshot_surface
127 
128  use mpi
129 
130  implicit none
131 
132  public :: init_snapshot_surface
133  public :: write_snapshot_surface
134  private :: write_snapshot_gmt_nat_fmt
135  private :: write_header_gmt_nat_fmt
136  private :: write_snapshot_vtk
137  public :: write_collection_vtk_surf
138  private :: collection_vtk_surface
139  public :: write_peak_ground_motion
140 
141  contains
142 
143 !
144 !
149 !***********************************************************************************************************************************************************************************
150  subroutine init_snapshot_surface()
151 !***********************************************************************************************************************************************************************************
152 
153  use mpi
154 
155  use mod_global_variables, only :&
156  tg_receiver_snapshot_quad&
157  ,ig_nquad_fsurf&
158  ,ig_quad_nnode&
159  ,ig_quadf_gll_glonum&
160  ,ig_quadf_gnode_glonum&
161  ,ig_nreceiver_snapshot&
162  ,ig_receiver_snapshot_glonum&
163  ,rg_receiver_snapshot_z&
164  ,ig_receiver_snapshot_locnum&
165  ,ig_receiver_snapshot_total_number&
166  ,ig_receiver_snapshot_mpi_shift&
167  ,ig_lst_unit&
168  ,ig_ngll&
169  ,ig_myrank&
170  ,error_stop&
171  ,ig_ncpu&
172  ,ig_receiver_snapshot_nx&
173  ,ig_receiver_snapshot_ny&
174  ,rg_mesh_xmin&
175  ,rg_mesh_xmax&
176  ,rg_mesh_ymin&
177  ,rg_mesh_ymax&
178  ,rg_receiver_snapshot_dxdy&
179  ,rg_receiver_snapshot_dx&
180  ,rg_receiver_snapshot_dy
181 
182  use mod_receiver , only :&
183  compute_info_quad_receiver&
184  ,is_inside_quad
185 
186  implicit none
187 
188  real , allocatable, dimension(:) :: rldum2
189  real , allocatable, dimension(:) :: rldum3
190  real :: rec_x
191  real :: rec_y
192  real :: rec_dmin
193 
194 
195  integer, allocatable, dimension(:) :: ildum1
196  integer :: nx
197  integer :: ny
198  integer :: ix
199  integer :: iy
200  integer :: iz
201 
202  integer :: rec_iel
203  integer :: rec_lgll
204  integer :: rec_mgll
205  integer :: irec
206  integer :: jrec
207  integer :: nrec
208 
209  integer :: icpu
210  integer :: iloc
211  integer :: ios
212 
213  logical, allocatable, dimension(:) :: is_rec_in_cpu
214  logical, allocatable, dimension(:) :: is_rec_in_all_cpu
215  logical :: is_inside
216  logical, dimension(ig_ncpu) :: is_inside_all_cpu
217 
218  character(len=255) :: info
219 
220  if (ig_myrank == 0) then
221  write(ig_lst_unit,'(" ",/,a)') "checking if receivers for snapshot are inside quadrangle elements..."
222  call flush(ig_lst_unit)
223  endif
224 
225 !
226 !---->adjust the space increment given in file *.cfg
227  rg_receiver_snapshot_dx = (rg_mesh_xmax - rg_mesh_xmin)/int((rg_mesh_xmax - rg_mesh_xmin)/rg_receiver_snapshot_dxdy)
228  rg_receiver_snapshot_dy = (rg_mesh_ymax - rg_mesh_ymin)/int((rg_mesh_ymax - rg_mesh_ymin)/rg_receiver_snapshot_dxdy)
229 
230 !
231 !---->set up the regular grid of receivers at the free surface
232  nx = int((rg_mesh_xmax - rg_mesh_xmin)/rg_receiver_snapshot_dx)+1
233  ny = int((rg_mesh_ymax - rg_mesh_ymin)/rg_receiver_snapshot_dy)+1
234 
235  ig_receiver_snapshot_nx = nx
236  ig_receiver_snapshot_ny = ny
237 
238  allocate(is_rec_in_cpu(nx*ny),stat=ios)
239  if (ios /= 0) then
240  write(info,'(a)') "error in subroutine init_snapshot_surface while allocating is_rec_in_cpu"
241  call error_stop(info)
242  else
243  do ix = 1,nx*ny
244  is_rec_in_cpu(ix) = .false.
245  enddo
246  endif
247 
248  allocate(is_rec_in_all_cpu(nx*ny),stat=ios)
249  if (ios /= 0) then
250  write(info,'(a)') "error in subroutine init_snapshot_surface while allocating is_rec_in_all_cpu"
251  call error_stop(info)
252  else
253  do ix = 1,nx*ny
254  is_rec_in_all_cpu(ix) = .false.
255  enddo
256  endif
257 
258  if (ig_myrank == 0) then
259  write(ig_lst_unit,'(2(a,F10.3))') " -->space increment dx,dy between receivers for snapshot = ",rg_receiver_snapshot_dx,",",rg_receiver_snapshot_dy
260  write(ig_lst_unit,'(3(a,i0))') " -->number of receivers for snapshot = ",nx,"*",ny," = ",nx * ny
261  call flush(ig_lst_unit)
262  endif
263 !
264 !---->first pass on nx*ny receivers to count the number of receivers that belong to cpu 'myrank'
265  irec = 0
266  nrec = 0
267 
268  do iy = 1,ny
269 
270  do ix = 1,nx
271 
272  irec = irec + 1
273  rec_x = rg_mesh_xmin + real(ix-1)*rg_receiver_snapshot_dx
274  rec_y = rg_mesh_ymax - real(iy-1)*rg_receiver_snapshot_dy
275 
276 !
277 !---------->find which quadrangle element contains the receiver irec
278  if (ig_nquad_fsurf > 0) then
279 
280  call is_inside_quad(rec_x,rec_y,ig_quadf_gnode_glonum,ig_nquad_fsurf,ig_quad_nnode,rec_dmin,rec_lgll,rec_mgll,rec_iel,is_inside)
281 
282  else
283 
284  is_inside = .false.
285 
286  endif
287 
288 !
289 !---------->cpu0 gather all is_inside to check: 1) if the receiver is inside a quadrangle element: if not abort computation 2) if the receiver is found by multiple cpus: if yes, affect the receiver to only one cpu
290  call mpi_gather(is_inside,1,mpi_logical,is_inside_all_cpu,1,mpi_logical,0,mpi_comm_world,ios)
291 
292  if (ig_myrank == 0) then
293 !
294 !------------->search the first cpu that contains the receiver
295  iloc = -1
296  do icpu = 1,ig_ncpu
297 
298  if (is_inside_all_cpu(icpu)) then
299  iloc = icpu
300  exit
301  endif
302 
303  enddo
304 
305 !
306 !------------->if the receiver is not inside any quadrangle element among all the cpus --> abort computation
307  if (iloc == -1) then
308 
309  write(info,'(a)') "error in subroutine init_snapshot_surface: receiver for snapshot not found over all cpus"
310  call error_stop(info)
311 
312  else
313 
314  do icpu = iloc+1,ig_ncpu
315  is_inside_all_cpu(icpu) = .false.
316  enddo
317 
318  endif
319 
320  endif
321 
322 !
323 !---------->cpu0 scatter array is_inside_all_cpu. The cpu that obtain the value is_inside = .true. will compute the receiver irec
324  call mpi_scatter(is_inside_all_cpu,1,mpi_logical,is_inside,1,mpi_logical,0,mpi_comm_world,ios)
325 
326  if (is_inside) then
327 
328  nrec = nrec + 1
329  is_rec_in_cpu(irec) = .true.
330 
331  endif
332 
333  if ( (ig_myrank == 0) .and. ( (irec == 1) .or. (mod(irec,5000) == 0) .or. (irec == nx*ny) ) ) then
334  write(ig_lst_unit,'(a,i10)') "checking receiver ",irec
335  call flush(ig_lst_unit)
336  endif
337 
338  enddo
339 
340  enddo
341 
342 !
343 !---->allocate array structure tg_receiver_snapshot_quad
344  ig_nreceiver_snapshot = nrec
345 
346  if (ig_nreceiver_snapshot > 0) then
347 
348  allocate(tg_receiver_snapshot_quad(ig_nreceiver_snapshot),stat=ios)
349 
350  if (ios /= 0) then
351 
352  write(info,'(a)') "error in subroutine init_snapshot_surface while allocating tg_receiver_snapshot_quad"
353  call error_stop(info)
354 
355  else
356 
357  do irec = 1,ig_nreceiver_snapshot
358 
359  tg_receiver_snapshot_quad(irec)%x = 0.0
360  tg_receiver_snapshot_quad(irec)%y = 0.0
361  tg_receiver_snapshot_quad(irec)%z = 0.0
362  tg_receiver_snapshot_quad(irec)%xi = 0.0
363  tg_receiver_snapshot_quad(irec)%eta = 0.0
364  tg_receiver_snapshot_quad(irec)%dmin = huge(rec_dmin)
365  tg_receiver_snapshot_quad(irec)%lag(:,:) = 0.0
366  tg_receiver_snapshot_quad(irec)%pgd_x = 0.0
367  tg_receiver_snapshot_quad(irec)%pgd_y = 0.0
368  tg_receiver_snapshot_quad(irec)%pgd_z = 0.0
369  tg_receiver_snapshot_quad(irec)%pgv_x = 0.0
370  tg_receiver_snapshot_quad(irec)%pgv_y = 0.0
371  tg_receiver_snapshot_quad(irec)%pgv_z = 0.0
372  tg_receiver_snapshot_quad(irec)%pgv_xyz = 0.0
373  tg_receiver_snapshot_quad(irec)%pga_x = 0.0
374  tg_receiver_snapshot_quad(irec)%pga_y = 0.0
375  tg_receiver_snapshot_quad(irec)%pga_z = 0.0
376  tg_receiver_snapshot_quad(irec)%gll(:,:) = 0
377  tg_receiver_snapshot_quad(irec)%cpu = 0
378  tg_receiver_snapshot_quad(irec)%iel = 0
379  tg_receiver_snapshot_quad(irec)%lgll = 0
380  tg_receiver_snapshot_quad(irec)%mgll = 0
381  tg_receiver_snapshot_quad(irec)%rglo = 0
382 
383  enddo
384 
385  endif
386 
387  endif
388 !
389 !---->second pass on receiver that belong to cpu 'myrank' and fill array structure tg_receiver_snapshot_quad
390  irec = 0
391  jrec = 0
392 
393  if ( (ig_nquad_fsurf > 0) .and. (ig_nreceiver_snapshot > 0) ) then
394 
395  do iy = 1,ny
396  do ix = 1,nx
397 
398  irec = irec + 1
399 
400  if (is_rec_in_cpu(irec)) then
401 
402  jrec = jrec + 1
403  rec_x = rg_mesh_xmin + real(ix-1)*rg_receiver_snapshot_dx
404  rec_y = rg_mesh_ymax - real(iy-1)*rg_receiver_snapshot_dy
405 
406 !
407 !---------------->find which quadrangle element contains the receiver irec
408  call is_inside_quad(rec_x,rec_y,ig_quadf_gnode_glonum,ig_nquad_fsurf,ig_quad_nnode,rec_dmin,rec_lgll,rec_mgll,rec_iel,is_inside)
409 
410 
411  if (is_inside) then
412 
413  tg_receiver_snapshot_quad(jrec)%x = rec_x
414  tg_receiver_snapshot_quad(jrec)%y = rec_y
415  tg_receiver_snapshot_quad(jrec)%dmin = rec_dmin
416  tg_receiver_snapshot_quad(jrec)%cpu = ig_myrank+1
417  tg_receiver_snapshot_quad(jrec)%iel = rec_iel
418  tg_receiver_snapshot_quad(jrec)%lgll = rec_lgll
419  tg_receiver_snapshot_quad(jrec)%mgll = rec_mgll
420  tg_receiver_snapshot_quad(jrec)%rglo = irec
421 
422  endif
423 
424  endif
425 
426  enddo
427  enddo
428 
429  endif
430 
431  deallocate(is_rec_in_cpu)
432 !
433 !---->compute receivers local coordinates and initialize PGx
434  if (ig_myrank == 0) then
435  write(ig_lst_unit,'(a)') " -->computing local coordinates of receivers"
436  call flush(ig_lst_unit)
437  endif
438 
439  do irec = 1,ig_nreceiver_snapshot
440 
441  call compute_info_quad_receiver(tg_receiver_snapshot_quad(irec))
442 
443  enddo
444 
445 !
446 !---->cpu0 gathers the number of receiver of all cpus. To avoid seg_fault at runtime with compilation option -check all, remove the if condition.
447  if (ig_myrank == 0) then
448  allocate(ig_receiver_snapshot_total_number(ig_ncpu))
449  allocate(ig_receiver_snapshot_mpi_shift(ig_ncpu))
450  endif
451 
452  call mpi_gather(ig_nreceiver_snapshot,1,mpi_integer,ig_receiver_snapshot_total_number,1,mpi_integer,0,mpi_comm_world,ios)
453 
454 !
455 !---->cpu0 prepares array ig_receiver_snapshot_mpi_shift for mpi_gatherv
456  if (ig_myrank == 0) then
457 
458  ig_receiver_snapshot_mpi_shift(1) = 0
459 
460  do icpu = 2,ig_ncpu
461  ig_receiver_snapshot_mpi_shift(icpu) = ig_receiver_snapshot_mpi_shift(icpu-1) + ig_receiver_snapshot_total_number(icpu-1)
462  enddo
463 
464  endif
465 !
466 !---->cpu0 gathers global number of all snapshot receivers of all cpus (because only cpu0 writes the *.grd or *.vts files)
467  allocate(ig_receiver_snapshot_glonum(nx*ny),stat=ios)
468  if (ios /= 0) then
469  write(info,'(a)') "error in subroutine init_snapshot_surface while allocating ig_receiver_snapshot_glonum"
470  call error_stop(info)
471  endif
472 !
473 !---->cpu0 gathers z-coordinate of all snapshot receivers of all cpus (because only cpu0 writes the *.grd or .*vts files)
474  allocate(rldum3(nx*ny),stat=ios)
475  if (ios /= 0) then
476  write(info,'(a)') "error in subroutine init_snapshot_surface while allocating rldum3"
477  call error_stop(info)
478  endif
479 !
480 !---->cpu0 gathers global number and z-coordinate of all snapshot receivers of all cpus
481  allocate(ildum1(ig_nreceiver_snapshot),stat=ios)
482  allocate(rldum2(ig_nreceiver_snapshot),stat=ios)
483 
484  if (ios /= 0) then
485 
486  write(info,'(a)') "error in subroutine init_snapshot_surface while allocating *dum*"
487  call error_stop(info)
488 
489  else
490 
491  do irec = 1,ig_nreceiver_snapshot
492  ildum1(irec) = tg_receiver_snapshot_quad(irec)%rglo
493  rldum2(irec) = tg_receiver_snapshot_quad(irec)%z
494  enddo
495 
496  endif
497 
498  call mpi_gatherv(ildum1 &
499  ,ig_nreceiver_snapshot &
500  ,mpi_integer &
501  ,ig_receiver_snapshot_glonum &
502  ,ig_receiver_snapshot_total_number &
503  ,ig_receiver_snapshot_mpi_shift &
504  ,mpi_integer &
505  ,0 &
506  ,mpi_comm_world &
507  ,ios )
508 
509  call mpi_gatherv(rldum2 &
510  ,ig_nreceiver_snapshot &
511  ,mpi_real &
512  ,rldum3 &
513  ,ig_receiver_snapshot_total_number &
514  ,ig_receiver_snapshot_mpi_shift &
515  ,mpi_real &
516  ,0 &
517  ,mpi_comm_world &
518  ,ios )
519 
520  deallocate(rldum2)
521 
522  if (ig_myrank == 0) then
523 
524  allocate(rg_receiver_snapshot_z(nx*ny),stat=ios)
525  if (ios /= 0) then
526  write(info,'(a)') "error in subroutine init_snapshot_surface while allocating rg_receiver_snapshot_z"
527  call error_stop(info)
528  endif
529 
530  allocate(ig_receiver_snapshot_locnum(nx*ny),stat=ios)
531  if (ios /= 0) then
532  write(info,'(a)') "error in subroutine init_snapshot_surface while allocating ig_receiver_snapshot_locnum"
533  call error_stop(info)
534  endif
535 
536  irec = 0
537 
538  do iy = 1,ny
539  do ix = 1,nx
540 
541  irec = irec + 1
542 
543  do iz = 1,nx*ny
544 
545  if (ig_receiver_snapshot_glonum(iz) == irec) then
546 
547  ig_receiver_snapshot_locnum(irec) = iz
548  rg_receiver_snapshot_z(irec) = rldum3(iz)
549 
550  exit
551 
552  endif
553 
554  enddo
555 
556  enddo
557  enddo
558 
559  endif
560 
561  if (ig_myrank == 0) then
562  write(ig_lst_unit,'(a)') "done"
563  call flush(ig_lst_unit)
564  endif
565 
566  return
567 
568 !***********************************************************************************************************************************************************************************
569  end subroutine init_snapshot_surface
570 !***********************************************************************************************************************************************************************************
571 
572 
573 !
574 !
578 !***********************************************************************************************************************************************************************************
579  subroutine write_snapshot_surface()
580 !***********************************************************************************************************************************************************************************
581 
582  use mpi
583 
584  use mod_global_variables, only :&
585  lg_snapshot_gmt&
586  ,lg_snapshot_vtk&
587  ,ig_ndof&
588  ,ig_ngll&
589  ,ig_idt&
590  ,cg_prefix&
591  ,lg_snapshot_displacement&
592  ,lg_snapshot_velocity&
593  ,lg_snapshot_acceleration&
594  ,ig_nreceiver_snapshot&
595  ,ig_snapshot_saving_incr&
596  ,tg_receiver_snapshot_quad&
597  ,rg_gll_displacement&
598  ,rg_gll_velocity&
599  ,rg_gll_acceleration&
600  ,ig_myrank
601 
602  use mod_gll_value , only : get_quad_gll_value
603 
604  use mod_lagrange , only : quad_lagrange_interp
605 
606  implicit none
607 
608  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: gll_dis
609  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: gll_vel
610  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: gll_acc
611 
612  real, dimension(ig_nreceiver_snapshot) :: ux
613  real, dimension(ig_nreceiver_snapshot) :: uy
614  real, dimension(ig_nreceiver_snapshot) :: uz
615  real, dimension(ig_nreceiver_snapshot) :: vx
616  real, dimension(ig_nreceiver_snapshot) :: vy
617  real, dimension(ig_nreceiver_snapshot) :: vz
618  real, dimension(ig_nreceiver_snapshot) :: ax
619  real, dimension(ig_nreceiver_snapshot) :: ay
620  real, dimension(ig_nreceiver_snapshot) :: az
621 
622  integer :: irec
623 
624  character(len=255) :: fname
625  character(len= 6) :: csnapshot
626  character(len= 4) :: vname
627 !
628 !
629 !*********************************************************************************************************************************
630 !---->before writing snapshots, displacement, velocity and acceleration are interpolated at the snapshot_receiver x,y,z location
631 !*********************************************************************************************************************************
632 
633  do irec = 1,ig_nreceiver_snapshot
634 
635 !
636 !------->get GLL nodes displacement, velocity and acceleration xyz-values from quadrangle element which contains the snapshot receiver
637  call get_quad_gll_value(rg_gll_displacement,tg_receiver_snapshot_quad(irec)%gll,gll_dis)
638  call get_quad_gll_value(rg_gll_velocity ,tg_receiver_snapshot_quad(irec)%gll,gll_vel)
639  call get_quad_gll_value(rg_gll_acceleration,tg_receiver_snapshot_quad(irec)%gll,gll_acc)
640 
641 !
642 !------->interpolate displacement at the snapshot receiver
643  call quad_lagrange_interp(gll_dis,tg_receiver_snapshot_quad(irec)%lag,ux(irec),uy(irec),uz(irec))
644 !
645 !------->compute PGD in x, y and z directions
646  if ( abs(ux(irec)) > tg_receiver_snapshot_quad(irec)%pgd_x ) tg_receiver_snapshot_quad(irec)%pgd_x = abs(ux(irec))
647  if ( abs(uy(irec)) > tg_receiver_snapshot_quad(irec)%pgd_y ) tg_receiver_snapshot_quad(irec)%pgd_y = abs(uy(irec))
648  if ( abs(uz(irec)) > tg_receiver_snapshot_quad(irec)%pgd_z ) tg_receiver_snapshot_quad(irec)%pgd_z = abs(uz(irec))
649 !
650 !------->interpolate velocity at the snapshot receiver
651  call quad_lagrange_interp(gll_vel,tg_receiver_snapshot_quad(irec)%lag,vx(irec),vy(irec),vz(irec))
652 !
653 !------->compute PGV in x, y and z directions
654  if ( abs(vx(irec)) > tg_receiver_snapshot_quad(irec)%pgv_x ) tg_receiver_snapshot_quad(irec)%pgv_x = abs(vx(irec))
655  if ( abs(vy(irec)) > tg_receiver_snapshot_quad(irec)%pgv_y ) tg_receiver_snapshot_quad(irec)%pgv_y = abs(vy(irec))
656  if ( abs(vz(irec)) > tg_receiver_snapshot_quad(irec)%pgv_z ) tg_receiver_snapshot_quad(irec)%pgv_z = abs(vz(irec))
657  if ( sqrt(vx(irec)**2 + vy(irec)**2 + vz(irec)**2) > tg_receiver_snapshot_quad(irec)%pgv_xyz ) tg_receiver_snapshot_quad(irec)%pgv_xyz = sqrt(vx(irec)**2 + vy(irec)**2 + vz(irec)**2)
658 !
659 !------->interpolate acceleration at the snapshot receiver
660  call quad_lagrange_interp(gll_acc,tg_receiver_snapshot_quad(irec)%lag,ax(irec),ay(irec),az(irec))
661 !
662 !------->compute PGA in x, y and z directions
663  if ( abs(ax(irec)) > tg_receiver_snapshot_quad(irec)%pga_x ) tg_receiver_snapshot_quad(irec)%pga_x = abs(ax(irec))
664  if ( abs(ay(irec)) > tg_receiver_snapshot_quad(irec)%pga_y ) tg_receiver_snapshot_quad(irec)%pga_y = abs(ay(irec))
665  if ( abs(az(irec)) > tg_receiver_snapshot_quad(irec)%pga_z ) tg_receiver_snapshot_quad(irec)%pga_z = abs(az(irec))
666 
667  enddo
668 !
669 !
670 !*********************************************************************************************************************************
671 !---->write snapshot in gmt or VTK format
672 !*********************************************************************************************************************************
673  write(csnapshot,'(i6.6)') ig_idt
674 !
675 !
676 !---->snapshot displacement
677  if ( lg_snapshot_displacement .and. (mod(ig_idt-1,ig_snapshot_saving_incr) == 0) ) then
678 
679  if (lg_snapshot_gmt) then
680 
681  fname = trim(cg_prefix)//".snapshot.ux."//trim(csnapshot)//".grd"
682  call write_snapshot_gmt_nat_fmt(fname,ux)
683 
684  fname = trim(cg_prefix)//".snapshot.uy."//trim(csnapshot)//".grd"
685  call write_snapshot_gmt_nat_fmt(fname,uy)
686 
687  fname = trim(cg_prefix)//".snapshot.uz."//trim(csnapshot)//".grd"
688  call write_snapshot_gmt_nat_fmt(fname,uz)
689 
690  endif
691 
692  if (lg_snapshot_vtk) then
693 
694  fname = trim(cg_prefix)//".snapshot.uxyz."//trim(csnapshot)//".vts"
695  vname = "disp"
696  call write_snapshot_vtk(fname,vname,ux,uy,uz)
697 
698  endif
699 
700  endif
701 !
702 !
703 !---->snapshot velocity
704  if ( lg_snapshot_velocity .and. (mod(ig_idt-1,ig_snapshot_saving_incr) == 0) ) then
705 
706  if (lg_snapshot_gmt) then
707 
708  fname = trim(cg_prefix)//".snapshot.vx."//trim(csnapshot)//".grd"
709  call write_snapshot_gmt_nat_fmt(fname,vx)
710 
711  fname = trim(cg_prefix)//".snapshot.vy."//trim(csnapshot)//".grd"
712  call write_snapshot_gmt_nat_fmt(fname,vy)
713 
714  fname = trim(cg_prefix)//".snapshot.vz."//trim(csnapshot)//".grd"
715  call write_snapshot_gmt_nat_fmt(fname,vz)
716 
717  endif
718 
719  if (lg_snapshot_vtk) then
720 
721  fname = trim(cg_prefix)//".snapshot.vxyz."//trim(csnapshot)//".vts"
722  vname = "velo"
723  call write_snapshot_vtk(fname,vname,vx,vy,vz)
724 
725  endif
726 
727  endif
728 !
729 !
730 !---->snapshot acceleration
731  if ( lg_snapshot_acceleration .and. (mod(ig_idt-1,ig_snapshot_saving_incr) == 0) ) then
732 
733  if (lg_snapshot_gmt) then
734 
735  fname = trim(cg_prefix)//".snapshot.ax."//trim(csnapshot)//".grd"
736  call write_snapshot_gmt_nat_fmt(fname,ax)
737 
738  fname = trim(cg_prefix)//".snapshot.ay."//trim(csnapshot)//".grd"
739  call write_snapshot_gmt_nat_fmt(fname,ay)
740 
741  fname = trim(cg_prefix)//".snapshot.az."//trim(csnapshot)//".grd"
742  call write_snapshot_gmt_nat_fmt(fname,az)
743 
744  endif
745 
746  if (lg_snapshot_vtk) then
747 
748  fname = trim(cg_prefix)//".snapshot.axyz."//trim(csnapshot)//".vts"
749  vname = "acce"
750  call write_snapshot_vtk(fname,vname,ax,ay,az)
751 
752  endif
753 
754  endif
755 
756  return
757 
758 !***********************************************************************************************************************************************************************************
759  end subroutine write_snapshot_surface
760 !***********************************************************************************************************************************************************************************
761 
762 !
763 !
768 !***********************************************************************************************************************************************************************************
769  subroutine write_snapshot_gmt_nat_fmt(fname,val)
770 !***********************************************************************************************************************************************************************************
771 
772  use mpi
773 
774  use mod_global_variables, only :&
775  ig_nreceiver_snapshot&
776  ,tg_receiver_snapshot_quad&
777  ,ig_receiver_snapshot_glonum&
778  ,ig_receiver_snapshot_total_number&
779  ,ig_receiver_snapshot_mpi_shift&
780  ,rg_receiver_snapshot_dx&
781  ,rg_receiver_snapshot_dy&
782  ,rg_mesh_xmax&
783  ,rg_mesh_ymax&
784  ,rg_mesh_xmin&
785  ,rg_mesh_ymin&
786  ,ig_receiver_snapshot_nx&
787  ,ig_receiver_snapshot_ny&
788  ,ig_myrank&
789  ,ig_ncpu
790 
791  use mod_init_memory, only : init_array_real
792 
793  implicit none
794 
795  character(len=255), intent(in) :: fname
796  real , intent(in) :: val(ig_nreceiver_snapshot)
797 
798  integer , parameter :: header_size = 892
799 
800  real :: val_max
801  real :: val_min
802  real :: val_max_all_cpu
803  real :: val_min_all_cpu
804  real, allocatable, dimension(:) :: val_all_cpu !allocated by cpu0 only. If compiled with -check all option, then this array must be allocated to all cpus to avoid seg_fault at runtime
805  real, allocatable, dimension(:) :: val_all_cpu_order
806 
807  integer :: nboctet_real
808  integer :: irec
809  integer :: rglo
810  integer :: desc
811  integer :: ios
812  integer(kind=mpi_offset_kind) :: pos
813  integer, dimension(mpi_status_size) :: statut
814 
815 !
816 !---->compute min and max. cpu 0 gather data for writing gmt header
817  if (ig_nreceiver_snapshot > 0) then
818  val_max = maxval(val)
819  val_min = minval(val)
820  else
821  val_min = 0.0
822  val_max = 0.0
823  endif
824  call mpi_reduce(val_max,val_max_all_cpu,1,mpi_real,mpi_max,0,mpi_comm_world,ios)
825  call mpi_reduce(val_min,val_min_all_cpu,1,mpi_real,mpi_min,0,mpi_comm_world,ios)
826 !
827 !---->open file 'fname' for snapshot
828  call mpi_file_open(mpi_comm_world,fname,mpi_mode_wronly + mpi_mode_create,mpi_info_null,desc,ios)
829  call mpi_type_size(mpi_real,nboctet_real,ios)
830 
831  if (ig_myrank == 0) then
832 !
833 !------>cpu 0 writes gmt grd file header
834  call write_header_gmt_nat_fmt(desc &
835  ,ig_receiver_snapshot_nx,ig_receiver_snapshot_ny &
836  ,rg_mesh_xmin,rg_mesh_xmax &
837  ,rg_mesh_ymin,rg_mesh_ymax &
838  ,val_min_all_cpu,val_max_all_cpu &
839  ,rg_receiver_snapshot_dx,rg_receiver_snapshot_dy )
840 
841  ios = init_array_real(val_all_cpu,ig_receiver_snapshot_nx*ig_receiver_snapshot_ny,"val_all_cpu")
842 
843  ios = init_array_real(val_all_cpu_order,ig_receiver_snapshot_nx*ig_receiver_snapshot_ny,"val_all_cpu_order")
844 
845  endif
846 
847 !
848 !---->cpu0 gathers array val to write gmt *.grd file
849  call mpi_gatherv(val &
850  ,ig_nreceiver_snapshot &
851  ,mpi_real &
852  ,val_all_cpu &
853  ,ig_receiver_snapshot_total_number &
854  ,ig_receiver_snapshot_mpi_shift &
855  ,mpi_real &
856  ,0 &
857  ,mpi_comm_world &
858  ,ios )
859 !
860 !---->cpu0 writes all receivers of all cpus using global numbering of snapshot receiver
861  if (ig_myrank == 0) then
862 
863 !
864 !------->order values gathered from all cpus so that cpu0 can write them at once
865  do irec = 1,ig_receiver_snapshot_nx*ig_receiver_snapshot_ny
866 
867  rglo = ig_receiver_snapshot_glonum(irec)
868 
869  val_all_cpu_order(rglo) = val_all_cpu(irec)
870 
871  enddo
872 
873 !
874 !------->cpu0 writes into GMT *.grd=bf file
875  pos = header_size
876  call mpi_file_write_at(desc,pos,val_all_cpu_order,ig_receiver_snapshot_nx*ig_receiver_snapshot_ny,mpi_real,statut,ios)
877 
878  endif
879 !
880 !---->close files
881  call mpi_file_close(desc,ios)
882 
883  return
884 
885 !***********************************************************************************************************************************************************************************
886  end subroutine write_snapshot_gmt_nat_fmt
887 !***********************************************************************************************************************************************************************************
888 
889 !
890 !
904 !***********************************************************************************************************************************************************************************
905  subroutine write_header_gmt_nat_fmt(desc,nx,ny,x_min,x_max,y_min,y_max,z_min,z_max,x_inc,y_inc)
906 !***********************************************************************************************************************************************************************************
907 
908  use mpi
909 
910  use mod_global_variables, only : ig_myrank
911 
912  implicit none
913 
914  integer, intent(in) :: desc
915  integer, intent(in) :: nx
916  integer, intent(in) :: ny
917  real , intent(in) :: x_min
918  real , intent(in) :: x_max
919  real , intent(in) :: y_min
920  real , intent(in) :: y_max
921  real , intent(in) :: z_min
922  real , intent(in) :: z_max
923  real , intent(in) :: x_inc
924  real , intent(in) :: y_inc
925 
926  integer, parameter :: grd_command_len = 320
927  integer, parameter :: grd_remark_len = 160
928  integer, parameter :: grd_title_len = 80
929  integer, parameter :: grd_unit_len = 80
930  integer, parameter :: grd_varname_len = 80
931  integer, parameter :: node_offset = 1 !pixel node registration if MNT file has an ESRI/ASCII format
932  real , parameter :: z_scale_factor = 1.0 !grd values must be multiplied by this
933  real , parameter :: z_add_offset = 0.0 !After scaling, add this
934 
935  double precision, dimension(10) :: ddum
936 
937  integer, dimension(3) :: idum
938  integer :: nboctet_double_precision
939  integer :: nboctet_char
940  integer :: nboctet_int
941  integer :: ios
942  integer(kind=mpi_offset_kind) :: pos
943  integer, dimension(mpi_status_size) :: statut
944 
945  character(len=GRD_UNIT_LEN) :: x_units !units in x-direction
946  character(len=GRD_UNIT_LEN) :: y_units !units in y-direction
947  character(len=GRD_UNIT_LEN) :: z_units !grid value units
948  character(len=GRD_TITLE_LEN) :: title !name of data set
949  character(len=GRD_COMMAND_LEN) :: command !name of generating command
950  character(len=GRD_REMARK_LEN) :: remark !comments on this data set
951 
952 
953  call mpi_type_size(mpi_double_precision,nboctet_double_precision,ios)
954  call mpi_type_size(mpi_character ,nboctet_char ,ios)
955  call mpi_type_size(mpi_integer ,nboctet_int ,ios)
956 
957  x_units = "m"//char(0)
958  y_units = "m"//char(0)
959  z_units = " "//char(0)
960  title = "EFISPEC3D snapshot"//char(0)
961  command = "mpi_file_write_at()"//char(0)
962  remark = "create by EFISPEC3D: http://efispec.free.fr"//char(0)
963 
964  idum(1) = nx
965  idum(2) = ny
966  idum(3) = node_offset
967 
968  ddum(1) = dble(x_min)
969  if (node_offset == 1) then
970  ddum(2) = dble(x_max)+dble(x_inc) !x_inc must be added to the size of the domain beause of pixel node registration (see variable NODE_OFFSET)
971  ddum(4) = dble(y_max)+dble(y_inc) !y_inc must be added to the size of the domain beause of pixel node registration (see variable NODE_OFFSET)
972  else
973  ddum(2) = dble(x_max)
974  ddum(4) = dble(y_max)
975  endif
976  ddum( 3) = dble(y_min)
977  ddum( 5) = dble(z_min)
978  ddum( 6) = dble(z_max)
979  ddum( 7) = dble(x_inc)
980  ddum( 8) = dble(y_inc)
981  ddum( 9) = dble(z_scale_factor)
982  ddum(10) = dble(z_add_offset)
983 
984  pos = 0
985  call mpi_file_write_at(desc,pos,idum,3,mpi_integer,statut,ios)
986 
987  pos = pos + 3*nboctet_int
988  call mpi_file_write_at(desc,pos,ddum,10,mpi_double_precision,statut,ios)
989 
990  pos = pos + 10*nboctet_double_precision
991  call mpi_file_write_at(desc,pos,x_units,grd_unit_len,mpi_character,statut,ios)
992 
993  pos = pos + nboctet_char*grd_unit_len
994  call mpi_file_write_at(desc,pos,y_units,grd_unit_len,mpi_character,statut,ios)
995 
996  pos = pos + nboctet_char*grd_unit_len
997  call mpi_file_write_at(desc,pos,z_units,grd_unit_len,mpi_character,statut,ios)
998 
999  pos = pos + nboctet_char*grd_unit_len
1000  call mpi_file_write_at(desc,pos,title,grd_title_len,mpi_character,statut,ios)
1001 
1002  pos = pos + nboctet_char*grd_title_len
1003  call mpi_file_write_at(desc,pos,command,grd_command_len,mpi_character,statut,ios)
1004 
1005  pos = pos + nboctet_char*grd_command_len
1006  call mpi_file_write_at(desc,pos,remark,grd_remark_len,mpi_character,statut,ios)
1007 
1008  return
1009 
1010 !***********************************************************************************************************************************************************************************
1011  end subroutine write_header_gmt_nat_fmt
1012 !***********************************************************************************************************************************************************************************
1013 
1014 !
1015 !
1023 !***********************************************************************************************************************************************************************************
1024  subroutine write_snapshot_vtk(fname,varname,valx,valy,valz)
1025 !***********************************************************************************************************************************************************************************
1026 
1027  use mpi
1028 
1029  use mod_global_variables, only :&
1030  ig_nreceiver_snapshot&
1031  ,ig_receiver_snapshot_total_number&
1032  ,ig_receiver_snapshot_mpi_shift&
1033  ,ig_receiver_snapshot_locnum&
1034  ,rg_receiver_snapshot_z&
1035  ,ig_myrank&
1036  ,error_stop&
1037  ,ig_receiver_snapshot_nx&
1038  ,ig_receiver_snapshot_ny&
1039  ,rg_receiver_snapshot_dx&
1040  ,rg_receiver_snapshot_dy&
1041  ,rg_mesh_xmin&
1042  ,rg_mesh_ymax
1043 
1044  use mod_vtk_io
1045 
1046  implicit none
1047 
1048  character(len=255), intent(in) :: fname
1049  character(len= 4), intent(in) :: varname
1050 
1051  real , intent(in) :: valx(ig_nreceiver_snapshot)
1052  real , intent(in) :: valy(ig_nreceiver_snapshot)
1053  real , intent(in) :: valz(ig_nreceiver_snapshot)
1054 
1055  real :: val_all_cpu(ig_receiver_snapshot_nx*ig_receiver_snapshot_ny)
1056  real :: val_all_cpu_order(ig_receiver_snapshot_nx*ig_receiver_snapshot_ny)
1057  real :: xrec(ig_receiver_snapshot_nx*ig_receiver_snapshot_ny)
1058  real :: yrec(ig_receiver_snapshot_nx*ig_receiver_snapshot_ny)
1059 
1060  integer :: ix
1061  integer :: iy
1062  integer :: irec
1063  integer :: irec_gatherv
1064  integer :: ios
1065 
1066  character(len=255) :: info
1067 
1068 !
1069 !---->cpu0 gathers array valx to be store in VTK file
1070  call mpi_gatherv(valx &
1071  ,ig_nreceiver_snapshot &
1072  ,mpi_real &
1073  ,val_all_cpu &
1074  ,ig_receiver_snapshot_total_number &
1075  ,ig_receiver_snapshot_mpi_shift &
1076  ,mpi_real &
1077  ,0 &
1078  ,mpi_comm_world &
1079  ,ios )
1080 !
1081 !---->cpu0 writes all receivers of all cpus using global numbering of snapshot receiver
1082  if (ig_myrank == 0) then
1083 
1084 !
1085 !------->compute x, y and z coordinates of receiver snapshot
1086  irec = 0
1087 
1088  do iy = 1,ig_receiver_snapshot_ny
1089  do ix = 1,ig_receiver_snapshot_nx
1090 
1091  irec = irec + 1
1092  irec_gatherv = ig_receiver_snapshot_locnum(irec)
1093 
1094  xrec(irec) = rg_mesh_xmin + real(ix-1)*rg_receiver_snapshot_dx
1095  yrec(irec) = rg_mesh_ymax - real(iy-1)*rg_receiver_snapshot_dy !y for VTK coordinate system
1096 
1097  val_all_cpu_order(irec) = val_all_cpu(irec_gatherv)
1098 
1099  enddo
1100  enddo
1101 
1102 !
1103 !------->initialize VTK file
1104  ios = vtk_ini_xml( &
1105  output_format = 'BINARY' &
1106  ,filename = trim(fname) &
1107  ,mesh_topology = 'StructuredGrid' &
1108  ,nx1 = 1 &
1109  ,nx2 = ig_receiver_snapshot_nx &
1110  ,ny1 = 1 &
1111  ,ny2 = ig_receiver_snapshot_ny &
1112  ,nz1 = 1 &
1113  ,nz2 = 1 )
1114 
1115  if (ios /= 0) then
1116  write(info,'(a)') "error in subroutine write_snapshot_vtk: VTK_INI_XML"
1117  call error_stop(info)
1118  endif
1119 
1120 !
1121 !------->store free surface geometry
1122  ios = vtk_geo_xml( &
1123  nx1 = 1 &
1124  ,nx2 = ig_receiver_snapshot_nx &
1125  ,ny1 = 1 &
1126  ,ny2 = ig_receiver_snapshot_ny &
1127  ,nz1 = 1 &
1128  ,nz2 = 1 &
1129  ,nn = ig_receiver_snapshot_nx*ig_receiver_snapshot_ny &
1130  ,x = xrec &
1131  ,y = yrec &
1132  ,z = rg_receiver_snapshot_z )
1133 
1134  if (ios /= 0) then
1135  write(info,'(a)') "error in subroutine write_snapshot_vtk: VTK_GEO_XML (init)"
1136  call error_stop(info)
1137  endif
1138 
1139 !
1140 !------->store x-component
1141  ios = vtk_dat_xml( &
1142  var_location = 'NODE'&
1143  ,var_block_action = 'OPEN')
1144 
1145  if (ios /= 0) then
1146  write(info,'(a)') "error in subroutine write_snapshot_vtk: VTK_DAT_XML (open)"
1147  call error_stop(info)
1148  endif
1149 
1150  ios = vtk_var_xml( &
1151  nc_nn = ig_receiver_snapshot_nx*ig_receiver_snapshot_ny &
1152  ,varname = "X-"//trim(varname) &
1153  ,var = val_all_cpu_order &
1154  )
1155 
1156  if (ios /= 0) then
1157  write(info,'(a)') "error in subroutine write_snapshot_vtk: VTK_VAR_XML X-component"
1158  call error_stop(info)
1159  endif
1160 
1161  endif
1162 
1163 !
1164 !---->cpu0 gathers array valy to be store in VTK file
1165  call mpi_gatherv(valy &
1166  ,ig_nreceiver_snapshot &
1167  ,mpi_real &
1168  ,val_all_cpu &
1169  ,ig_receiver_snapshot_total_number &
1170  ,ig_receiver_snapshot_mpi_shift &
1171  ,mpi_real &
1172  ,0 &
1173  ,mpi_comm_world &
1174  ,ios )
1175 
1176 
1177  if (ig_myrank == 0) then
1178 
1179 !
1180 !------->order data from gatherv
1181  irec = 0
1182  do iy = 1,ig_receiver_snapshot_ny
1183  do ix = 1,ig_receiver_snapshot_nx
1184 
1185  irec = irec + 1
1186  irec_gatherv = ig_receiver_snapshot_locnum(irec)
1187  val_all_cpu_order(irec) = val_all_cpu(irec_gatherv)
1188 
1189  enddo
1190  enddo
1191 
1192 !
1193 !------->store y-component
1194  ios = vtk_var_xml( &
1195  nc_nn = ig_receiver_snapshot_nx*ig_receiver_snapshot_ny &
1196  ,varname = "Y-"//trim(varname) &
1197  ,var = val_all_cpu_order &
1198  )
1199  if (ios /= 0) then
1200  write(info,'(a)') "error in subroutine write_snapshot_vtk: VTK_VAR_XML Y-component"
1201  call error_stop(info)
1202  endif
1203 
1204  endif
1205 
1206 !
1207 !---->cpu0 gathers array valz to be store in VTK file
1208  call mpi_gatherv(valz &
1209  ,ig_nreceiver_snapshot &
1210  ,mpi_real &
1211  ,val_all_cpu &
1212  ,ig_receiver_snapshot_total_number &
1213  ,ig_receiver_snapshot_mpi_shift &
1214  ,mpi_real &
1215  ,0 &
1216  ,mpi_comm_world &
1217  ,ios )
1218 
1219  if (ig_myrank == 0) then
1220 
1221 !
1222 !------->order data from gatherv
1223  irec = 0
1224  do iy = 1,ig_receiver_snapshot_ny
1225  do ix = 1,ig_receiver_snapshot_nx
1226 
1227  irec = irec + 1
1228  irec_gatherv = ig_receiver_snapshot_locnum(irec)
1229  val_all_cpu_order(irec) = val_all_cpu(irec_gatherv)
1230 
1231  enddo
1232  enddo
1233 !
1234 !------->store z-component
1235  ios = vtk_var_xml( &
1236  nc_nn = ig_receiver_snapshot_nx*ig_receiver_snapshot_ny &
1237  ,varname = "Z-"//trim(varname) &
1238  ,var = val_all_cpu_order &
1239  )
1240 
1241  if (ios /= 0) then
1242  write(info,'(a)') "error in subroutine write_snapshot_vtk: VTK_VAR_XML Z-component"
1243  call error_stop(info)
1244  endif
1245 
1246 !
1247 !------->finalize VTK file
1248  ios = vtk_dat_xml( &
1249  var_location = 'NODE' &
1250  ,var_block_action = 'CLOSE')
1251 
1252  if (ios /= 0) then
1253  write(info,'(a)') "error in subroutine write_snapshot_vtk: VTK_DAT_XML (close)"
1254  call error_stop(info)
1255  endif
1256 
1257  ios = vtk_geo_xml()
1258 
1259  if (ios /= 0) then
1260  write(info,'(a)') "error in subroutine write_snapshot_vtk: VTK_GEO_XML (finalize)"
1261  call error_stop(info)
1262  endif
1263 
1264  ios = vtk_end_xml()
1265 
1266  if (ios /= 0) then
1267  write(info,'(a)') "error in subroutine write_snapshot_vtk: VTK_END_XML"
1268  call error_stop(info)
1269  endif
1270 
1271  endif
1272 
1273  return
1274 !***********************************************************************************************************************************************************************************
1275  end subroutine write_snapshot_vtk
1276 !***********************************************************************************************************************************************************************************
1277 
1278 !
1279 !
1283 !***********************************************************************************************************************************************************************************
1284  subroutine write_collection_vtk_surf()
1285 !***********************************************************************************************************************************************************************************
1286 
1287  use mod_global_variables, only :&
1288  lg_snapshot_displacement&
1289  ,lg_snapshot_velocity&
1290  ,lg_snapshot_acceleration&
1291  ,lg_snapshot_volume_displacement&
1292  ,lg_snapshot_volume_velocity&
1293  ,lg_snapshot_volume_acceleration
1294 
1295  implicit none
1296 
1297  if (lg_snapshot_displacement) then
1298  call collection_vtk_surface("surface.snapshot.uxyz")
1299  endif
1300 
1301  if (lg_snapshot_velocity) then
1302  call collection_vtk_surface("surface.snapshot.vxyz")
1303  endif
1304 
1305  if (lg_snapshot_acceleration) then
1306  call collection_vtk_surface("surface.snapshot.axyz")
1307  endif
1308 
1309  return
1310 !***********************************************************************************************************************************************************************************
1311  end subroutine write_collection_vtk_surf
1312 !***********************************************************************************************************************************************************************************
1313 
1314 !
1315 !
1319 !***********************************************************************************************************************************************************************************
1320  subroutine collection_vtk_surface(varname)
1321 !***********************************************************************************************************************************************************************************
1322 
1323  use mod_global_variables, only :&
1324  get_newunit&
1325  ,cg_prefix&
1326  ,rg_dt&
1327  ,ig_ndt&
1328  ,ig_snapshot_saving_incr
1329 
1330  implicit none
1331 
1332  character(len=*), intent(in) :: varname
1333 
1334  real :: time
1335 
1336  integer :: istep
1337  integer :: myunit
1338 
1339  character(len=255) :: fname
1340  character(len=6 ) :: csnapshot
1341 
1342  open(unit=get_newunit(myunit),file=trim(cg_prefix)//".collection."//trim(varname)//".pvd")
1343 
1344  write(unit=myunit,fmt='(a)') "<?xml version=""1.0""?>"
1345  write(unit=myunit,fmt='(a)') "<VTKFile type=""Collection"" version=""0.1"" byte_order=""BigEndian"">"
1346  write(unit=myunit,fmt='(a)') " <Collection>"
1347 
1348  do istep = 1,ig_ndt,ig_snapshot_saving_incr
1349 
1350  write(csnapshot,'(I6.6)') istep
1351 
1352  time = real(istep-1)*rg_dt
1353 
1354  fname = trim(cg_prefix)//"."//trim(varname)//"."//trim(csnapshot)//".vts"
1355 
1356  write(unit=myunit,fmt='(a,E14.7,3a)') " <DataSet timestep=""",time,""" group="""" part=""0"" file=""",trim(fname),"""/>"
1357 
1358  enddo
1359 
1360  write(unit=myunit,fmt='(a)') " </Collection>"
1361  write(unit=myunit,fmt='(a)') "</VTKFile>"
1362 
1363  close(myunit)
1364 
1365  return
1366 !***********************************************************************************************************************************************************************************
1367  end subroutine collection_vtk_surface
1368 !***********************************************************************************************************************************************************************************
1369 
1370 !
1371 !
1374 !***********************************************************************************************************************************************************************************
1375  subroutine write_peak_ground_motion()
1376 !***********************************************************************************************************************************************************************************
1377 
1378  use mpi
1379 
1380  use mod_global_variables, only :&
1381  tg_receiver_snapshot_quad&
1382  ,ig_myrank&
1383  ,ig_nreceiver_snapshot&
1384  ,cg_prefix&
1385  ,ig_lst_unit
1386 
1387  implicit none
1388 
1389  real, dimension(ig_nreceiver_snapshot) :: pgd_x
1390  real, dimension(ig_nreceiver_snapshot) :: pgd_y
1391  real, dimension(ig_nreceiver_snapshot) :: pgd_z
1392  real, dimension(ig_nreceiver_snapshot) :: pgv_x
1393  real, dimension(ig_nreceiver_snapshot) :: pgv_y
1394  real, dimension(ig_nreceiver_snapshot) :: pgv_z
1395  real, dimension(ig_nreceiver_snapshot) :: pga_x
1396  real, dimension(ig_nreceiver_snapshot) :: pga_y
1397  real, dimension(ig_nreceiver_snapshot) :: pga_z
1398 
1399  integer :: irec
1400 
1401  character(len=255) :: fname
1402 
1403  if (ig_myrank == 0) then
1404  write(ig_lst_unit,'(" ",/,a)') "writing snapshots of peak ground motion..."
1405  call flush(ig_lst_unit)
1406  endif
1407 
1408  do irec = 1,ig_nreceiver_snapshot
1409  pgd_x(irec) = tg_receiver_snapshot_quad(irec)%pgd_x
1410  enddo
1411 
1412  do irec = 1,ig_nreceiver_snapshot
1413  pgd_y(irec) = tg_receiver_snapshot_quad(irec)%pgd_y
1414  enddo
1415 
1416  do irec = 1,ig_nreceiver_snapshot
1417  pgd_z(irec) = tg_receiver_snapshot_quad(irec)%pgd_z
1418  enddo
1419 
1420  do irec = 1,ig_nreceiver_snapshot
1421  pgv_x(irec) = tg_receiver_snapshot_quad(irec)%pgv_x
1422  enddo
1423 
1424  do irec = 1,ig_nreceiver_snapshot
1425  pgv_y(irec) = tg_receiver_snapshot_quad(irec)%pgv_y
1426  enddo
1427 
1428  do irec = 1,ig_nreceiver_snapshot
1429  pgv_z(irec) = tg_receiver_snapshot_quad(irec)%pgv_z
1430  enddo
1431 
1432  do irec = 1,ig_nreceiver_snapshot
1433  pga_x(irec) = tg_receiver_snapshot_quad(irec)%pga_x
1434  enddo
1435 
1436  do irec = 1,ig_nreceiver_snapshot
1437  pga_y(irec) = tg_receiver_snapshot_quad(irec)%pga_y
1438  enddo
1439 
1440  do irec = 1,ig_nreceiver_snapshot
1441  pga_z(irec) = tg_receiver_snapshot_quad(irec)%pga_z
1442  enddo
1443 
1444 !
1445 !---->write peak ground motion in GMT format
1446  fname = trim(cg_prefix)//".snapshot.PGDx.grd"
1447  call write_snapshot_gmt_nat_fmt(fname,pgd_x)
1448 
1449  fname = trim(cg_prefix)//".snapshot.PGDy.grd"
1450  call write_snapshot_gmt_nat_fmt(fname,pgd_y)
1451 
1452  fname = trim(cg_prefix)//".snapshot.PGDz.grd"
1453  call write_snapshot_gmt_nat_fmt(fname,pgd_z)
1454 
1455  fname = trim(cg_prefix)//".snapshot.PGVx.grd"
1456  call write_snapshot_gmt_nat_fmt(fname,pgv_x)
1457 
1458  fname = trim(cg_prefix)//".snapshot.PGVy.grd"
1459  call write_snapshot_gmt_nat_fmt(fname,pgv_y)
1460 
1461  fname = trim(cg_prefix)//".snapshot.PGVz.grd"
1462  call write_snapshot_gmt_nat_fmt(fname,pgv_z)
1463 
1464  fname = trim(cg_prefix)//".snapshot.PGAx.grd"
1465  call write_snapshot_gmt_nat_fmt(fname,pga_x)
1466 
1467  fname = trim(cg_prefix)//".snapshot.PGAy.grd"
1468  call write_snapshot_gmt_nat_fmt(fname,pga_y)
1469 
1470  fname = trim(cg_prefix)//".snapshot.PGAz.grd"
1471  call write_snapshot_gmt_nat_fmt(fname,pga_z)
1472 
1473 !
1474 !---->write peak ground motion in VTK format
1475  fname = trim(cg_prefix)//".snapshot.PGDxyz.vts"
1476  call write_snapshot_vtk(fname,"displacement",pgd_x,pgd_y,pgd_z)
1477 
1478  fname = trim(cg_prefix)//".snapshot.PGVxyz.vts"
1479  call write_snapshot_vtk(fname,"velocity" ,pgv_x,pgv_y,pgv_z)
1480 
1481  fname = trim(cg_prefix)//".snapshot.PGAxyz.vts"
1482  call write_snapshot_vtk(fname,"acceleration",pga_x,pga_y,pga_z)
1483 
1484  if (ig_myrank == 0) then
1485  write(ig_lst_unit,'(a)') "done"
1486  call flush(ig_lst_unit)
1487  endif
1488 
1489  return
1490 
1491 !***********************************************************************************************************************************************************************************
1492  end subroutine write_peak_ground_motion
1493 !***********************************************************************************************************************************************************************************
1494 
1495 end module mod_snapshot_surface
mod_init_memory
This module contains subroutines to allocate arrays and to compute an approximation of the total RAM ...
Definition: module_init_memory.f90:126
mod_receiver
This module contains subroutines to compute information about receivers and to write receivers&#39; time ...
Definition: module_receiver.f90:126
mod_snapshot_surface::write_snapshot_vtk
subroutine, private write_snapshot_vtk(fname, varname, valx, valy, valz)
This subroutine writes structured grid snapshot in binary VTK xml format.
Definition: module_snapshot_surface.f90:1024
mod_vtk_io::vtk_end_xml
integer(i4p) function, public vtk_end_xml()
The VTK_END_XML function finalizes opened files.
Definition: module_vtk_io.f90:960
mod_vtk_io::vtk_ini_xml
integer(i4p) function, public vtk_ini_xml(output_format, filename, mesh_topology, nx1, nx2, ny1, ny2, nz1, nz2)
The VTK_INI_XML function is used for initializing file. This function must be the first to be called...
Definition: module_vtk_io.f90:288
mod_vtk_io::vtk_geo_xml
overloading of VTK_GEO_XML
Definition: module_vtk_io.f90:157
mod_snapshot_surface::write_snapshot_gmt_nat_fmt
subroutine, private write_snapshot_gmt_nat_fmt(fname, val)
This subroutine writes structured grid snapshot in binary native GMT format.
Definition: module_snapshot_surface.f90:769
mod_global_variables
This module defines all global variables of EFISPEC3D. Scalar variables are initialized directly in t...
Definition: module_global_variables.f90:128
mod_init_memory::init_array_real
Interface init_array_real to redirect allocation to n-rank arrays.
Definition: module_init_memory.f90:148
mod_snapshot_surface::init_snapshot_surface
subroutine, public init_snapshot_surface()
This subroutine generates a structured grid of receivers on the free surface.
Definition: module_snapshot_surface.f90:150
mod_receiver::compute_info_quad_receiver
subroutine, public compute_info_quad_receiver(tlxinf)
This subroutine computes local coordinates of a receiver in the quadrangle element it belongs...
Definition: module_receiver.f90:917
mod_snapshot_surface::write_header_gmt_nat_fmt
subroutine, private write_header_gmt_nat_fmt(desc, nx, ny, x_min, x_max, y_min, y_max, z_min, z_max, x_inc, y_inc)
This subroutine writes header of binary native GMT files.
Definition: module_snapshot_surface.f90:905
mod_vtk_io::vtk_dat_xml
integer(i4p) function, public vtk_dat_xml(var_location, var_block_action)
The VTK_DAT_XML function opens or closes CellData/PointData structures.
Definition: module_vtk_io.f90:695
mod_snapshot_surface::collection_vtk_surface
subroutine, private collection_vtk_surface(varname)
This subroutine write Paraview collection file *.pvd of VTK xml snapshot files *.vts.
Definition: module_snapshot_surface.f90:1320
mod_vtk_io
This module programmed by S. Zaghi (https://github.com/szaghi/Lib_VTK_IO) contains functions to write...
Definition: module_vtk_io.f90:130
mod_lagrange::quad_lagrange_interp
subroutine, public quad_lagrange_interp(gll_values, lag, x, y, z)
This subroutine interpolates GLL nodes -values of a quadrangle element using Lagrange polynomials...
Definition: module_lagrange.f90:365
mod_snapshot_surface::write_peak_ground_motion
subroutine, public write_peak_ground_motion()
This subroutine writes peak ground motions files in GMT or VTK formats.
Definition: module_snapshot_surface.f90:1375
mod_snapshot_surface::write_collection_vtk_surf
subroutine, public write_collection_vtk_surf()
This subroutine selects which ParaView collection files should be written (displacement, velocity and/or acceleration collections) depending on value of variables mod_global_variables::lg_snapshot_displacement, mod_global_variables::lg_snapshot_velocity and mod_global_variables::lg_snapshot_acceleration.
Definition: module_snapshot_surface.f90:1284
mod_gll_value
This module contains subroutines to get GLL nodes values from global GLL nodes numbering.
Definition: module_gll_value.f90:126
mod_gll_value::get_quad_gll_value
subroutine, public get_quad_gll_value(gll_all_values, gll_global_number, gll_selected_values)
subroutine to get quadrangle element GLL values from an array ordered by global GLL nodes numbering...
Definition: module_gll_value.f90:200
mod_receiver::is_inside_quad
subroutine, public is_inside_quad(x, y, quad_gnode, nquad, quad_nnode, dmin, lgll, mgll, quad_num, is_inside)
This subroutine finds which quadrangle element contains a receiver.
Definition: module_receiver.f90:1143
mod_lagrange
This module contains functions to compute Lagrange polynomials and its derivatives; and to interpolat...
Definition: module_lagrange.f90:126
mod_vtk_io::vtk_var_xml
overloading of VTK_VAR_XML
Definition: module_vtk_io.f90:166
mod_snapshot_surface
This module contains subroutines to compute and write snapshots of the free surface (either in GMT or...
Definition: module_snapshot_surface.f90:126
mod_snapshot_surface::write_snapshot_surface
subroutine, public write_snapshot_surface()
This subroutine computes and writes -displacements, velocities and accelerations at receivers used fo...
Definition: module_snapshot_surface.f90:579