From 2935eef43b739cb0f02548b88b992409537cc4f2 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Mon, 19 Apr 2021 23:21:21 +0300 Subject: [PATCH] renames xstore_dummy/.. to xstore_unique/.. for generate databases --- .../create_mass_matrices.f90 | 6 +- .../create_regions_mesh.f90 | 80 ++-- .../fault_generate_databases.f90 | 66 ++-- src/generate_databases/generate_databases.f90 | 26 +- .../generate_databases_par.F90 | 20 +- src/generate_databases/get_MPI.f90 | 6 +- .../get_absorbing_boundary.F90 | 64 +-- .../get_coupling_surfaces.f90 | 44 +-- src/generate_databases/get_model.F90 | 18 +- src/generate_databases/model_1d_cascadia.f90 | 2 +- src/generate_databases/model_1d_prem.f90 | 2 +- .../model_external_values.F90 | 24 +- src/generate_databases/save_arrays_solver.F90 | 50 +-- .../save_arrays_solver_adios.F90 | 31 +- src/generate_databases/setup_color_perm.f90 | 6 +- src/generate_databases/setup_mesh.f90 | 23 +- .../input_output/Teleseismic_IO.f90 | 6 +- .../input_output/input_output_mod.f90 | 12 +- .../projection_on_FD_grid_mod.f90 | 10 +- .../regularization/regularization_SEM_mod.f90 | 53 ++- src/shared/get_element_face.f90 | 42 +- src/shared/get_jacobian_boundaries.f90 | 374 +++++++++--------- src/shared/write_VTK_data.f90 | 60 +-- src/specfem3D/setup_sources_receivers.f90 | 4 + src/tomography/save_external_bin_m_up.f90 | 30 +- .../Paraview/create_slice_VTK.f90 | 22 +- 26 files changed, 551 insertions(+), 530 deletions(-) diff --git a/src/generate_databases/create_mass_matrices.f90 b/src/generate_databases/create_mass_matrices.f90 index 3b9570f2c..e4ddf40ad 100644 --- a/src/generate_databases/create_mass_matrices.f90 +++ b/src/generate_databases/create_mass_matrices.f90 @@ -251,8 +251,8 @@ subroutine create_mass_matrices_ocean_load(nglob,nspec,ibool) if (TOPOGRAPHY) then ! takes elevation from topography file - xloc = xstore_dummy(iglobnum) - yloc = ystore_dummy(iglobnum) + xloc = xstore_unique(iglobnum) + yloc = ystore_unique(iglobnum) call get_topo_bathy_elevation(xloc,yloc,loc_elevation, & itopo_bathy,NX_TOPO,NY_TOPO) @@ -262,7 +262,7 @@ subroutine create_mass_matrices_ocean_load(nglob,nspec,ibool) else ! takes elevation from z-coordinate of mesh point - elevation = zstore_dummy(iglobnum) + elevation = zstore_unique(iglobnum) endif diff --git a/src/generate_databases/create_regions_mesh.f90 b/src/generate_databases/create_regions_mesh.f90 index 08f1e3a69..ee8d3449f 100644 --- a/src/generate_databases/create_regions_mesh.f90 +++ b/src/generate_databases/create_regions_mesh.f90 @@ -131,12 +131,12 @@ subroutine create_regions_mesh() nodes_coords_ext_mesh,nnodes_ext_mesh, & elmnts_ext_mesh,nelmnts_ext_mesh) endif - ! at this point (xyz)store_dummy are still open + ! at this point (xyz)store_unique are still open if (.not. PARALLEL_FAULT) then call fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, & xstore,ystore,zstore,nspec,nglob) endif - ! this closes (xyz)store_dummy + ! this closes (xyz)store_unique endif @@ -147,15 +147,15 @@ subroutine create_regions_mesh() write(IMAIN,*) ' ...preparing MPI interfaces ' call flush_IMAIN() endif - call get_MPI_interface(nglob_dummy,nspec,ibool) + call get_MPI_interface(nglob_unique,nspec,ibool) ! setting up parallel fault if (PARALLEL_FAULT .and. ANY_FAULT) then call synchronize_all() - !at this point (xyz)store_dummy are still open + !at this point (xyz)store_unique are still open call fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, & xstore,ystore,zstore,nspec,nglob) - ! this closes (xyz)store_dummy + ! this closes (xyz)store_unique endif ! sets up absorbing/free surface boundaries @@ -247,7 +247,7 @@ subroutine create_regions_mesh() write(IMAIN,*) ' ...creating C-PML damping profiles ' call flush_IMAIN() endif - call pml_set_local_dampingcoeff(xstore_dummy,ystore_dummy,zstore_dummy) + call pml_set_local_dampingcoeff(xstore_unique,ystore_unique,zstore_unique) endif ! creates mass matrix @@ -257,7 +257,7 @@ subroutine create_regions_mesh() write(IMAIN,*) ' ...creating mass matrix ' call flush_IMAIN() endif - call create_mass_matrices(nglob_dummy,nspec,ibool,PML_CONDITIONS,STACEY_ABSORBING_CONDITIONS) + call create_mass_matrices(nglob_unique,nspec,ibool,PML_CONDITIONS,STACEY_ABSORBING_CONDITIONS) ! saves the binary mesh files call synchronize_all() @@ -301,7 +301,7 @@ subroutine create_regions_mesh() call synchronize_all() ! computes the approximate amount of memory needed to run the solver - call memory_eval(nspec,nglob_dummy,maxval(nibool_interfaces_ext_mesh),num_interfaces_ext_mesh, & + call memory_eval(nspec,nglob_unique,maxval(nibool_interfaces_ext_mesh),num_interfaces_ext_mesh, & APPROXIMATE_OCEAN_LOAD,memory_size) call max_all_dp(memory_size, max_memory_size) @@ -313,8 +313,8 @@ subroutine create_regions_mesh() write(IMAIN,*) ' ...checking mesh resolution' call flush_IMAIN() endif - call check_mesh_resolution(nspec,nglob_dummy, & - ibool,xstore_dummy,ystore_dummy,zstore_dummy, & + call check_mesh_resolution(nspec,nglob_unique, & + ibool,xstore_unique,ystore_unique,zstore_unique, & ispec_is_acoustic,ispec_is_elastic,ispec_is_poroelastic, & kappastore,mustore,rhostore, & phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI, & @@ -346,7 +346,7 @@ subroutine create_regions_mesh() deallocate(rho_vpI,rho_vpII,rho_vsI) deallocate(rhoarraystore,kappaarraystore,etastore,phistore,tortstore,permstore) - if (.not. SAVE_MOHO_MESH) deallocate(xstore_dummy,ystore_dummy,zstore_dummy) + if (.not. SAVE_MOHO_MESH) deallocate(xstore_unique,ystore_unique,zstore_unique) if (ACOUSTIC_SIMULATION) deallocate(rmass_acoustic) if (ELASTIC_SIMULATION) deallocate(rmass) @@ -1086,28 +1086,30 @@ subroutine crm_ext_setup_indexing(ibool, & endif ! unique global point locations - nglob_dummy = nglob - allocate(xstore_dummy(nglob_dummy),stat=ier) + nglob_unique = nglob ! note: nglob => NGLOB_AB + + ! coordinates for global points + allocate(xstore_unique(nglob_unique),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 813') - xstore_dummy(:) = 0._CUSTOM_REAL + xstore_unique(:) = 0._CUSTOM_REAL - allocate(ystore_dummy(nglob_dummy),stat=ier) + allocate(ystore_unique(nglob_unique),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 814') - ystore_dummy(:) = 0._CUSTOM_REAL + ystore_unique(:) = 0._CUSTOM_REAL - allocate(zstore_dummy(nglob_dummy),stat=ier) + allocate(zstore_unique(nglob_unique),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 815') if (ier /= 0) stop 'error in allocate' - zstore_dummy(:) = 0._CUSTOM_REAL + zstore_unique(:) = 0._CUSTOM_REAL do ispec = 1, nspec do k = 1, NGLLZ do j = 1, NGLLY do i = 1, NGLLX iglobnum = ibool(i,j,k,ispec) - xstore_dummy(iglobnum) = real(xstore(i,j,k,ispec),kind=CUSTOM_REAL) - ystore_dummy(iglobnum) = real(ystore(i,j,k,ispec),kind=CUSTOM_REAL) - zstore_dummy(iglobnum) = real(zstore(i,j,k,ispec),kind=CUSTOM_REAL) + xstore_unique(iglobnum) = real(xstore(i,j,k,ispec),kind=CUSTOM_REAL) + ystore_unique(iglobnum) = real(ystore(i,j,k,ispec),kind=CUSTOM_REAL) + zstore_unique(iglobnum) = real(zstore(i,j,k,ispec),kind=CUSTOM_REAL) enddo enddo enddo @@ -1177,9 +1179,9 @@ subroutine crm_setup_moho(nspec,nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, & reshape( (/ 1,2,2, NGLLX,2,2, 2,1,2, 2,NGLLY,2, 2,2,1, 2,2,NGLLZ /),(/3,6/)) ! top ! temporary arrays for passing information - allocate(iglob_is_surface(nglob_dummy),stat=ier) + allocate(iglob_is_surface(nglob_unique),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 816') - allocate(iglob_normals(NDIM,nglob_dummy),stat=ier) + allocate(iglob_normals(NDIM,nglob_unique),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 817') if (ier /= 0) stop 'error allocating array iglob_is_surface' @@ -1204,8 +1206,8 @@ subroutine crm_setup_moho(nspec,nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, & ! sets face id of reference element associated with this face call get_element_face_id(ispec,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & iface) ! ijk indices of GLL points for face id @@ -1213,7 +1215,7 @@ subroutine crm_setup_moho(nspec,nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, & ! weighted jacobian and normal call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) @@ -1223,8 +1225,8 @@ subroutine crm_setup_moho(nspec,nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, & do j = 1,NGLLY do i = 1,NGLLX call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & normal_face(:,i,j)) enddo enddo @@ -1290,9 +1292,9 @@ subroutine crm_setup_moho(nspec,nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, & counter = counter+1 ! reference corner coordinates - xcoord(icorner) = xstore_dummy(iglob) - ycoord(icorner) = ystore_dummy(iglob) - zcoord(icorner) = zstore_dummy(iglob) + xcoord(icorner) = xstore_unique(iglob) + ycoord(icorner) = ystore_unique(iglob) + zcoord(icorner) = zstore_unique(iglob) endif enddo @@ -1305,7 +1307,7 @@ subroutine crm_setup_moho(nspec,nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, & ! re-computes face infos ! weighted jacobian and normal call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) @@ -1315,8 +1317,8 @@ subroutine crm_setup_moho(nspec,nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, & do j = 1,NGLLZ do i = 1,NGLLX call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & normal_face(:,i,j) ) enddo enddo @@ -1330,8 +1332,8 @@ subroutine crm_setup_moho(nspec,nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, & ! determines whether normal points into element or not (top/bottom distinction) call get_element_face_normal_idirect(ispec,iface,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & normal,idirect) ! takes moho surface element id given by id on midpoint @@ -1501,7 +1503,7 @@ subroutine crm_setup_inner_outer_elemnts(nspec,ibool,SAVE_MESH_FILES) if (ier /= 0) stop 'error allocating array ispec_is_inner' ! temporary array - allocate(iglob_is_inner(nglob_dummy),stat=ier) + allocate(iglob_is_inner(nglob_unique),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 827') if (ier /= 0) stop 'error allocating temporary array iglob_is_inner' @@ -1532,8 +1534,8 @@ subroutine crm_setup_inner_outer_elemnts(nspec,ibool,SAVE_MESH_FILES) if (SAVE_MESH_FILES .and. DEBUG) then filename = prname(1:len_trim(prname))//'ispec_is_inner' - call write_VTK_data_elem_l(nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + call write_VTK_data_elem_l(nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,ibool, & ispec_is_inner,filename) endif diff --git a/src/generate_databases/fault_generate_databases.f90 b/src/generate_databases/fault_generate_databases.f90 index 59229aa13..03c91dd08 100644 --- a/src/generate_databases/fault_generate_databases.f90 +++ b/src/generate_databases/fault_generate_databases.f90 @@ -213,7 +213,7 @@ subroutine fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, & ! checks if anything to do on this fault if (fault_db(iflt)%nspec == 0) cycle - !NOTE: the small fault opening in *_dummy does not affect this subroutine (see get_element_face_id) + !NOTE: the small fault opening in *_unique does not affect this subroutine (see get_element_face_id) call setup_iface(fault_db(iflt),nnodes_ext_mesh,nodes_coords_ext_mesh,nspec,nglob,ibool) ! saving GLL indices for each fault face, needed for ibulks @@ -227,7 +227,7 @@ subroutine fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, & ! from global indices on the fault to global indices on the bulk call setup_ibulks(fault_db(iflt),ibool,nspec) - ! close the fault in (xyz)store_dummy + ! close the fault in (xyz)store_unique call close_fault(fault_db(iflt)) call setup_Kelvin_Voigt_eta(fault_db(iflt),nspec) @@ -248,7 +248,7 @@ end subroutine fault_setup subroutine setup_iface(fdb,nnodes_ext_mesh,nodes_coords_ext_mesh,nspec,nglob,ibool) - use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy + use create_regions_mesh_ext_par, only: xstore_unique,ystore_unique,zstore_unique implicit none type(fault_db_type), intent(inout) :: fdb @@ -273,7 +273,7 @@ subroutine setup_iface(fdb,nnodes_ext_mesh,nodes_coords_ext_mesh,nspec,nglob,ibo enddo call get_element_face_id(fdb%ispec1(e),xcoord,ycoord,zcoord, & ibool,nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy, & + xstore_unique,ystore_unique,zstore_unique, & fdb%iface1(e)) ! side 2 do icorner = 1,NGNOD2D @@ -283,7 +283,7 @@ subroutine setup_iface(fdb,nnodes_ext_mesh,nodes_coords_ext_mesh,nspec,nglob,ibo enddo call get_element_face_id(fdb%ispec2(e),xcoord,ycoord,zcoord, & ibool,nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy, & + xstore_unique,ystore_unique,zstore_unique, & fdb%iface2(e)) enddo @@ -479,12 +479,12 @@ subroutine setup_ibulks(fdb,ibool,nspec) end subroutine setup_ibulks !============================================================================================================= -! We only close *store_dummy. *store is close already in create_regions_mesh.f90 -! Fortunately only *store_dummy is needed to compute jacobians and normals +! We only close *store_unique. *store is close already in create_regions_mesh.f90 +! Fortunately only *store_unique is needed to compute jacobians and normals subroutine close_fault(fdb) - use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy + use create_regions_mesh_ext_par, only: xstore_unique,ystore_unique,zstore_unique implicit none type(fault_db_type), intent(inout) :: fdb @@ -497,21 +497,21 @@ subroutine close_fault(fdb) K1 = fdb%ibulk1(i) K2 = fdb%ibulk2(i) - x1 = dble(xstore_dummy(K1)) - x2 = dble(xstore_dummy(K2)) + x1 = dble(xstore_unique(K1)) + x2 = dble(xstore_unique(K2)) - y1 = dble(ystore_dummy(K1)) - y2 = dble(ystore_dummy(K2)) + y1 = dble(ystore_unique(K1)) + y2 = dble(ystore_unique(K2)) - z1 = dble(zstore_dummy(K1)) - z2 = dble(zstore_dummy(K2)) + z1 = dble(zstore_unique(K1)) + z2 = dble(zstore_unique(K2)) - xstore_dummy(K1) = real(0.5d0 *(x1 + x2),kind=CUSTOM_REAL) - xstore_dummy(K2) = xstore_dummy(K1) - ystore_dummy(K1) = real(0.5d0 *(y1 + y2),kind=CUSTOM_REAL) - ystore_dummy(K2) = ystore_dummy(K1) - zstore_dummy(K1) = real(0.5d0 *(z1 + z2),kind=CUSTOM_REAL) - zstore_dummy(K2) = zstore_dummy(K1) + xstore_unique(K1) = real(0.5d0 *(x1 + x2),kind=CUSTOM_REAL) + xstore_unique(K2) = xstore_unique(K1) + ystore_unique(K1) = real(0.5d0 *(y1 + y2),kind=CUSTOM_REAL) + ystore_unique(K2) = ystore_unique(K1) + zstore_unique(K1) = real(0.5d0 *(z1 + z2),kind=CUSTOM_REAL) + zstore_unique(K2) = zstore_unique(K1) enddo end subroutine close_fault @@ -520,7 +520,7 @@ end subroutine close_fault subroutine save_fault_xyzcoord_ibulk(fdb) - use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy + use create_regions_mesh_ext_par, only: xstore_unique,ystore_unique,zstore_unique implicit none type(fault_db_type), intent(inout) :: fdb @@ -546,13 +546,13 @@ subroutine save_fault_xyzcoord_ibulk(fdb) do i = 1, fdb%nglob K1 = fdb%ibulk1(i) K2 = fdb%ibulk2(i) - fdb%xcoordbulk1(i) = xstore_dummy(K1) - fdb%ycoordbulk1(i) = ystore_dummy(K1) - fdb%zcoordbulk1(i) = zstore_dummy(K1) + fdb%xcoordbulk1(i) = xstore_unique(K1) + fdb%ycoordbulk1(i) = ystore_unique(K1) + fdb%zcoordbulk1(i) = zstore_unique(K1) - fdb%xcoordbulk2(i) = xstore_dummy(K2) - fdb%ycoordbulk2(i) = ystore_dummy(K2) - fdb%zcoordbulk2(i) = zstore_dummy(K2) + fdb%xcoordbulk2(i) = xstore_unique(K2) + fdb%ycoordbulk2(i) = ystore_unique(K2) + fdb%zcoordbulk2(i) = zstore_unique(K2) enddo end subroutine save_fault_xyzcoord_ibulk @@ -562,7 +562,7 @@ end subroutine save_fault_xyzcoord_ibulk subroutine setup_normal_jacobian(fdb,ibool,nspec,nglob) - use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy, & + use create_regions_mesh_ext_par, only: xstore_unique,ystore_unique,zstore_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz @@ -602,14 +602,14 @@ subroutine setup_normal_jacobian(fdb,ibool,nspec,nglob) iglob_corners_ref(icorner) = ibool(i,j,k,ispec) ! reference corner coordinates - xcoord(icorner) = xstore_dummy(iglob_corners_ref(icorner)) - ycoord(icorner) = ystore_dummy(iglob_corners_ref(icorner)) - zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner)) + xcoord(icorner) = xstore_unique(iglob_corners_ref(icorner)) + ycoord(icorner) = ystore_unique(iglob_corners_ref(icorner)) + zcoord(icorner) = zstore_unique(iglob_corners_ref(icorner)) enddo ! gets face GLL 2Djacobian, weighted from element face call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) @@ -620,7 +620,7 @@ subroutine setup_normal_jacobian(fdb,ibool,nspec,nglob) ! directs normals such that they point outwards of element call get_element_face_normal(ispec,iface_ref,xcoord,ycoord,zcoord, & ibool,nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy, & + xstore_unique,ystore_unique,zstore_unique, & normal_face(:,i,j) ) enddo enddo diff --git a/src/generate_databases/generate_databases.f90 b/src/generate_databases/generate_databases.f90 index f026df3d8..24bf9109e 100644 --- a/src/generate_databases/generate_databases.f90 +++ b/src/generate_databases/generate_databases.f90 @@ -195,21 +195,21 @@ program xgenerate_databases ! timing double precision, external :: wtime -! MPI initialization + ! MPI initialization call init_mpi() -! sizeprocs returns number of processes started (should be equal to NPROC). -! myrank is the rank of each process, between 0 and NPROC-1. -! as usual in MPI, process 0 is in charge of coordinating everything -! and also takes care of the main output + ! sizeprocs returns number of processes started (should be equal to NPROC). + ! myrank is the rank of each process, between 0 and NPROC-1. + ! as usual in MPI, process 0 is in charge of coordinating everything + ! and also takes care of the main output call world_size(sizeprocs) call world_rank(myrank) -! open main output file, only written to by process 0 + ! open main output file, only written to by process 0 if (myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) & open(unit=IMAIN,file=trim(OUTPUT_FILES)//'/output_generate_databases.txt',status='unknown') -! get MPI starting time + ! get MPI starting time time_start = wtime() if (myrank == 0) then @@ -225,10 +225,10 @@ program xgenerate_databases call flush_IMAIN() endif -! read the parameter file + ! read the parameter file call read_parameters() -! reads topography and bathymetry file + ! reads topography and bathymetry file call read_topography() if (myrank == 0) then @@ -245,24 +245,24 @@ program xgenerate_databases call adios_setup() endif -! reads Databases files + ! reads Databases files if (ADIOS_FOR_DATABASES) then call read_partition_files_adios() else call read_partition_files() endif -! external mesh creation + ! external mesh creation call setup_mesh() -! finalize mesher + ! finalize mesher call finalize_databases() if (ADIOS_ENABLED) then call adios_cleanup() endif -! MPI finish + ! MPI finish call finalize_mpi() end program xgenerate_databases diff --git a/src/generate_databases/generate_databases_par.F90 b/src/generate_databases/generate_databases_par.F90 index 685fa4335..6c62610ab 100644 --- a/src/generate_databases/generate_databases_par.F90 +++ b/src/generate_databases/generate_databases_par.F90 @@ -72,7 +72,10 @@ module generate_databases_par double precision :: max_memory_size,max_memory_size_request ! this for all the regions - integer :: NSPEC_AB,NGLOB_AB + ! number of elements + integer :: NSPEC_AB + ! number of (unique) global nodes + integer :: NGLOB_AB integer :: NSPEC2D_BOTTOM,NSPEC2D_TOP @@ -169,19 +172,20 @@ module create_regions_mesh_ext_par implicit none ! global point coordinates - real(kind=CUSTOM_REAL), dimension(:), allocatable :: xstore_dummy - real(kind=CUSTOM_REAL), dimension(:), allocatable :: ystore_dummy - real(kind=CUSTOM_REAL), dimension(:), allocatable :: zstore_dummy - integer :: nglob_dummy + real(kind=CUSTOM_REAL), dimension(:), allocatable :: xstore_unique + real(kind=CUSTOM_REAL), dimension(:), allocatable :: ystore_unique + real(kind=CUSTOM_REAL), dimension(:), allocatable :: zstore_unique + ! number of unique global points + integer :: nglob_unique ! ****** ! added by Ping Tong (TP / Tong Ping) for the FK3D calculation ! not used so far... ! material properties at global points for the elastic case - !real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhostore_dummy - !real(kind=CUSTOM_REAL), dimension(:), allocatable :: kappastore_dummy - !real(kind=CUSTOM_REAL), dimension(:), allocatable :: mustore_dummy + !real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhostore_unique + !real(kind=CUSTOM_REAL), dimension(:), allocatable :: kappastore_unique + !real(kind=CUSTOM_REAL), dimension(:), allocatable :: mustore_unique ! ****** ! Gauss-Lobatto-Legendre points and weights of integration diff --git a/src/generate_databases/get_MPI.f90 b/src/generate_databases/get_MPI.f90 index fff344463..4f2329eb1 100644 --- a/src/generate_databases/get_MPI.f90 +++ b/src/generate_databases/get_MPI.f90 @@ -151,9 +151,9 @@ subroutine get_MPI_interface(nglob,nspec,ibool) do ilocnum = 1, num_interface_points iglob = ibool_interfaces_ext_mesh(ilocnum,iinterface) ! we will use double precision locations to find/sort MPI points - xp(ilocnum) = dble(xstore_dummy(iglob)) - yp(ilocnum) = dble(ystore_dummy(iglob)) - zp(ilocnum) = dble(zstore_dummy(iglob)) + xp(ilocnum) = dble(xstore_unique(iglob)) + yp(ilocnum) = dble(ystore_unique(iglob)) + zp(ilocnum) = dble(zstore_unique(iglob)) enddo ! sorts (lexicographically?) ibool_interfaces_ext_mesh and updates value diff --git a/src/generate_databases/get_absorbing_boundary.F90 b/src/generate_databases/get_absorbing_boundary.F90 index 01b7231af..f125f9b35 100644 --- a/src/generate_databases/get_absorbing_boundary.F90 +++ b/src/generate_databases/get_absorbing_boundary.F90 @@ -149,8 +149,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! sets face id of reference element associated with this face call get_element_face_id(ispec,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy,iface) + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,iface) if (COUPLE_WITH_INJECTION_TECHNIQUE .or. MESH_A_CHUNK_OF_THE_EARTH) then iboun(iface,ispec) = .true. @@ -161,7 +161,7 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! weighted jacobian and normal call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) @@ -172,15 +172,15 @@ subroutine get_absorbing_boundary(nspec,ibool, & do i = 1,NGLLX lnormal(:) = normal_face(:,i,j) call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & lnormal ) normal_face(:,i,j) = lnormal(:) ! writes out xmin boundary point locations if ( COUPLE_WITH_INJECTION_TECHNIQUE .and. INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_DSM ) then - write(123,'(i10,3f20.10)') ispec, xstore_dummy(ibool(i,j,1,ispec)), ystore_dummy(ibool(i,j,1,ispec)), & - zstore_dummy(ibool(i,j,1,ispec)) + write(123,'(i10,3f20.10)') ispec, xstore_unique(ibool(i,j,1,ispec)), ystore_unique(ibool(i,j,1,ispec)), & + zstore_unique(ibool(i,j,1,ispec)) endif enddo enddo @@ -236,8 +236,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! sets face id of reference element associated with this face call get_element_face_id(ispec,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy,iface) + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,iface) if (COUPLE_WITH_INJECTION_TECHNIQUE .or. MESH_A_CHUNK_OF_THE_EARTH) then iboun(iface,ispec) = .true. @@ -248,7 +248,7 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! weighted jacobian and normal call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) @@ -259,8 +259,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & do i = 1,NGLLX lnormal(:) = normal_face(:,i,j) call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & lnormal ) normal_face(:,i,j) = lnormal(:) enddo @@ -313,8 +313,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! sets face id of reference element associated with this face call get_element_face_id(ispec,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy,iface) + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,iface) if (COUPLE_WITH_INJECTION_TECHNIQUE .or. MESH_A_CHUNK_OF_THE_EARTH) then iboun(iface,ispec) = .true. @@ -325,7 +325,7 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! weighted jacobian and normal call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) @@ -336,8 +336,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & do i = 1,NGLLY lnormal(:) = normal_face(:,i,j) call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & lnormal ) normal_face(:,i,j) = lnormal(:) enddo @@ -390,8 +390,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! sets face id of reference element associated with this face call get_element_face_id(ispec,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy,iface) + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,iface) if (COUPLE_WITH_INJECTION_TECHNIQUE .or. MESH_A_CHUNK_OF_THE_EARTH) then iboun(iface,ispec) = .true. @@ -402,7 +402,7 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! weighted jacobian and normal call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) @@ -413,8 +413,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & do i = 1,NGLLY lnormal(:) = normal_face(:,i,j) call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & lnormal ) normal_face(:,i,j) = lnormal(:) enddo @@ -467,8 +467,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! sets face id of reference element associated with this face call get_element_face_id(ispec,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy,iface) + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,iface) if (COUPLE_WITH_INJECTION_TECHNIQUE .or. MESH_A_CHUNK_OF_THE_EARTH) then iboun(iface,ispec) = .true. @@ -479,7 +479,7 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! weighted jacobian and normal call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) @@ -490,8 +490,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & do i = 1,NGLLX lnormal(:) = normal_face(:,i,j) call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & lnormal ) normal_face(:,i,j) = lnormal(:) enddo @@ -570,8 +570,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! sets face id of reference element associated with this face call get_element_face_id(ispec,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy,iface) + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,iface) if (COUPLE_WITH_INJECTION_TECHNIQUE .or. MESH_A_CHUNK_OF_THE_EARTH) then iboun(iface,ispec) = .true. @@ -582,7 +582,7 @@ subroutine get_absorbing_boundary(nspec,ibool, & ! weighted jacobian and normal call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) @@ -593,8 +593,8 @@ subroutine get_absorbing_boundary(nspec,ibool, & do i = 1,NGLLX lnormal(:) = normal_face(:,i,j) call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & lnormal ) normal_face(:,i,j) = lnormal(:) enddo diff --git a/src/generate_databases/get_coupling_surfaces.f90 b/src/generate_databases/get_coupling_surfaces.f90 index 6bc94b555..4e2866c12 100644 --- a/src/generate_databases/get_coupling_surfaces.f90 +++ b/src/generate_databases/get_coupling_surfaces.f90 @@ -63,13 +63,13 @@ subroutine get_coupling_surfaces(nspec,ibool) num_coupling_el_po_faces = 0 ! sets flags for acoustic / elastic / poroelastic on global points - allocate(acoustic_flag(nglob_dummy),stat=ier) + allocate(acoustic_flag(nglob_unique),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 690') if (ier /= 0) stop 'error allocating array acoustic_flag' - allocate(elastic_flag(nglob_dummy),stat=ier) + allocate(elastic_flag(nglob_unique),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 691') if (ier /= 0) stop 'error allocating array elastic_flag' - allocate(poroelastic_flag(nglob_dummy),stat=ier) + allocate(poroelastic_flag(nglob_unique),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 692') if (ier /= 0) stop 'error allocating array poroelastic_flag' @@ -132,7 +132,7 @@ subroutine get_coupling_surfaces(nspec,ibool) ! sums acoustic flags if (ACOUSTIC_SIMULATION) then - call assemble_MPI_scalar_i_blocking(NPROC,nglob_dummy,acoustic_flag, & + call assemble_MPI_scalar_i_blocking(NPROC,nglob_unique,acoustic_flag, & num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, & nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy, & my_neighbors_ext_mesh) @@ -140,7 +140,7 @@ subroutine get_coupling_surfaces(nspec,ibool) ! sums elastic flags if (ELASTIC_SIMULATION) then - call assemble_MPI_scalar_i_blocking(NPROC,nglob_dummy,elastic_flag, & + call assemble_MPI_scalar_i_blocking(NPROC,nglob_unique,elastic_flag, & num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, & nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy, & my_neighbors_ext_mesh) @@ -148,7 +148,7 @@ subroutine get_coupling_surfaces(nspec,ibool) ! sums poroelastic flags if (POROELASTIC_SIMULATION) then - call assemble_MPI_scalar_i_blocking(NPROC,nglob_dummy,poroelastic_flag, & + call assemble_MPI_scalar_i_blocking(NPROC,nglob_unique,poroelastic_flag, & num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, & nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh_dummy, & my_neighbors_ext_mesh) @@ -210,7 +210,7 @@ subroutine get_coupling_surfaces_ac_el(nspec,ibool,elastic_flag) ! arrays with the mesh integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool - integer,dimension(nglob_dummy),intent(in) :: elastic_flag + integer,dimension(nglob_unique),intent(in) :: elastic_flag ! local parameters ! (assumes NGLLX=NGLLY=NGLLZ) @@ -267,7 +267,7 @@ subroutine get_coupling_surfaces_ac_el(nspec,ibool,elastic_flag) tmp_normal(:,:,:) = 0.0 tmp_jacobian2Dw(:,:) = 0.0 - allocate(mask_ibool(nglob_dummy),stat=ier) + allocate(mask_ibool(nglob_unique),stat=ier) if (ier /= 0) call exit_MPI(myrank,'error allocating array 698') if (ier /= 0) stop 'error allocating array mask_ibool' mask_ibool(:) = .false. @@ -291,7 +291,7 @@ subroutine get_coupling_surfaces_ac_el(nspec,ibool,elastic_flag) ! takes indices of corners of reference face call get_element_corners(ispec,iface_ref,xcoord,ycoord,zcoord,iglob_corners_ref, & - ibool,nspec,nglob_dummy,xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique,xstore_unique,ystore_unique,zstore_unique, & iface_all_corner_ijk) ! checks if face is has an elastic side @@ -318,7 +318,7 @@ subroutine get_coupling_surfaces_ac_el(nspec,ibool,elastic_flag) ! gets face GLL 2Djacobian, weighted from element face call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) @@ -329,8 +329,8 @@ subroutine get_coupling_surfaces_ac_el(nspec,ibool,elastic_flag) do i=1,NGLLX ! directs normals such that they point outwards of element call get_element_face_normal(ispec,iface_ref,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & normal_face(:,i,j) ) ! makes sure that it always points away from acoustic element, ! otherwise switch direction @@ -425,7 +425,7 @@ subroutine get_coupling_surfaces_ac_poro(nspec,ibool,acoustic_flag) ! arrays with the mesh integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool - integer,dimension(nglob_dummy),intent(in) :: acoustic_flag + integer,dimension(nglob_unique),intent(in) :: acoustic_flag ! local parameters ! (assumes NGLLX=NGLLY=NGLLZ) @@ -490,7 +490,7 @@ subroutine get_coupling_surfaces_ac_poro(nspec,ibool,acoustic_flag) ! takes indices of corners of reference face call get_element_corners(ispec,iface_ref,xcoord,ycoord,zcoord,iglob_corners_ref, & - ibool,nspec,nglob_dummy,xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique,xstore_unique,ystore_unique,zstore_unique, & iface_all_corner_ijk) ! checks if face has acoustic side @@ -504,7 +504,7 @@ subroutine get_coupling_surfaces_ac_poro(nspec,ibool,acoustic_flag) ! gets face GLL 2Djacobian, weighted from element face call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) @@ -515,8 +515,8 @@ subroutine get_coupling_surfaces_ac_poro(nspec,ibool,acoustic_flag) do i = 1,NGLLX ! directs normals such that they point outwards of element call get_element_face_normal(ispec,iface_ref,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & normal_face(:,i,j) ) ! reverse the sign, we know we are in a poroelastic element ! thus, pointing outwards of acoustic element @@ -594,7 +594,7 @@ subroutine get_coupling_surfaces_el_poro(nspec,ibool,elastic_flag) ! arrays with the mesh integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool - integer,dimension(nglob_dummy),intent(in) :: elastic_flag + integer,dimension(nglob_unique),intent(in) :: elastic_flag ! local parameters ! (assumes NGLLX=NGLLY=NGLLZ) @@ -675,7 +675,7 @@ subroutine get_coupling_surfaces_el_poro(nspec,ibool,elastic_flag) ! takes indices of corners of reference face call get_element_corners(ispec,iface_ref,xcoord,ycoord,zcoord,iglob_corners_ref, & - ibool,nspec,nglob_dummy,xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique,xstore_unique,ystore_unique,zstore_unique, & iface_all_corner_ijk) ! checks if face has elastic side @@ -689,7 +689,7 @@ subroutine get_coupling_surfaces_el_poro(nspec,ibool,elastic_flag) ! gets face GLL 2Djacobian, weighted from element face call get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + xstore_unique,ystore_unique,zstore_unique,ibool,nglob_unique, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) @@ -700,8 +700,8 @@ subroutine get_coupling_surfaces_el_poro(nspec,ibool,elastic_flag) do i = 1,NGLLX ! directs normals such that they point outwards of poroelastic element call get_element_face_normal(ispec,iface_ref,xcoord,ycoord,zcoord, & - ibool,nspec,nglob_dummy, & - xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & normal_face(:,i,j) ) enddo enddo diff --git a/src/generate_databases/get_model.F90 b/src/generate_databases/get_model.F90 index 507b98da6..4cde1820b 100644 --- a/src/generate_databases/get_model.F90 +++ b/src/generate_databases/get_model.F90 @@ -123,9 +123,9 @@ subroutine get_model() if ((NGLLX == 5) .and. (NGLLY == 5) .and. (NGLLZ == 5)) then ! gets xyz coordinates of GLL point iglob = ibool(3,3,3,1) - xmesh = xstore_dummy(iglob) - ymesh = ystore_dummy(iglob) - zmesh = zstore_dummy(iglob) + xmesh = xstore_unique(iglob) + ymesh = ystore_unique(iglob) + zmesh = zstore_unique(iglob) call model_coupled_FindLayer(xmesh,ymesh,zmesh) else stop 'bad number of GLL points for coupling with an external code' @@ -143,9 +143,9 @@ subroutine get_model() if (COUPLE_WITH_INJECTION_TECHNIQUE .or. MESH_A_CHUNK_OF_THE_EARTH) then if (INJECTION_TECHNIQUE_TYPE /= INJECTION_TECHNIQUE_IS_FK) then iglob = ibool(3,3,3,ispec) - xmesh = xstore_dummy(iglob) - ymesh = ystore_dummy(iglob) - zmesh = zstore_dummy(iglob) + xmesh = xstore_unique(iglob) + ymesh = ystore_unique(iglob) + zmesh = zstore_unique(iglob) call model_coupled_FindLayer(xmesh,ymesh,zmesh) endif endif @@ -203,9 +203,9 @@ subroutine get_model() ! gets xyz coordinates of GLL point iglob = ibool(i,j,k,ispec) - xmesh = xstore_dummy(iglob) - ymesh = ystore_dummy(iglob) - zmesh = zstore_dummy(iglob) + xmesh = xstore_unique(iglob) + ymesh = ystore_unique(iglob) + zmesh = zstore_unique(iglob) !! VM VM for coupling with DSM !! find the layer in which the middle of the element is located diff --git a/src/generate_databases/model_1d_cascadia.f90 b/src/generate_databases/model_1d_cascadia.f90 index fa759cbd7..79c00b676 100644 --- a/src/generate_databases/model_1d_cascadia.f90 +++ b/src/generate_databases/model_1d_cascadia.f90 @@ -60,7 +60,7 @@ subroutine model_1D_cascadia(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten,qkappa_atten) distmin = HUGEVAL elevation = 0.0 call get_topo_elevation_free_closest(x,y,elevation,distmin, & - nspec,nglob_dummy,ibool,xstore_dummy,ystore_dummy,zstore_dummy, & + nspec,nglob_unique,ibool,xstore_unique,ystore_unique,zstore_unique, & num_free_surface_faces,free_surface_ispec,free_surface_ijk) ! depth in Z-direction diff --git a/src/generate_databases/model_1d_prem.f90 b/src/generate_databases/model_1d_prem.f90 index 6bc359c7e..e58591510 100644 --- a/src/generate_databases/model_1d_prem.f90 +++ b/src/generate_databases/model_1d_prem.f90 @@ -83,7 +83,7 @@ subroutine model_1D_prem_iso(xmesh,ymesh,zmesh,rho_prem,vp_prem,vs_prem,qmu_atte distmin = HUGEVAL elevation = 0.0 call get_topo_elevation_free_closest(xloc,yloc,elevation,distmin, & - nspec,nglob_dummy,ibool,xstore_dummy,ystore_dummy,zstore_dummy, & + nspec,nglob_unique,ibool,xstore_unique,ystore_unique,zstore_unique, & num_free_surface_faces,free_surface_ispec,free_surface_ijk) ! depth in Z-direction diff --git a/src/generate_databases/model_external_values.F90 b/src/generate_databases/model_external_values.F90 index da11e54b3..39c84ecfc 100644 --- a/src/generate_databases/model_external_values.F90 +++ b/src/generate_databases/model_external_values.F90 @@ -132,8 +132,8 @@ subroutine model_external_values(xmesh,ymesh,zmesh,ispec,rho,vp,vs,qkappa_atten, use generate_databases_par, only: nspec => NSPEC_AB,ibool,HUGEVAL,TINYVAL,IDOMAIN_ELASTIC,CUSTOM_REAL - use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy, & - num_free_surface_faces,free_surface_ijk,free_surface_ispec,nglob_dummy + use create_regions_mesh_ext_par, only: xstore_unique,ystore_unique,zstore_unique, & + num_free_surface_faces,free_surface_ijk,free_surface_ispec,nglob_unique implicit none @@ -177,26 +177,26 @@ subroutine model_external_values(xmesh,ymesh,zmesh,ispec,rho,vp,vs,qkappa_atten, ! gets element midpoint location (if needed) idummy = ispec !iglob = ibool(MIDX,MIDY,MIDZ,ispec) - !mid_x = xstore_dummy(iglob) - !mid_y = ystore_dummy(iglob) - !mid_z = zstore_dummy(iglob) + !mid_x = xstore_unique(iglob) + !mid_y = ystore_unique(iglob) + !mid_z = zstore_unique(iglob) ! note: z coordinate will be negative below surface ! convention is z-axis points up ! model dimensions - xmin = 0._CUSTOM_REAL ! minval(xstore_dummy) - xmax = 134000._CUSTOM_REAL ! maxval(xstore_dummy) - ymin = 0._CUSTOM_REAL !minval(ystore_dummy) - ymax = 134000._CUSTOM_REAL ! maxval(ystore_dummy) - zmin = 0._CUSTOM_REAL ! minval(zstore_dummy) - zmax = 60000._CUSTOM_REAL ! maxval(zstore_dummy) + xmin = 0._CUSTOM_REAL ! minval(xstore_unique) + xmax = 134000._CUSTOM_REAL ! maxval(xstore_unique) + ymin = 0._CUSTOM_REAL !minval(ystore_unique) + ymax = 134000._CUSTOM_REAL ! maxval(ystore_unique) + zmin = 0._CUSTOM_REAL ! minval(zstore_unique) + zmax = 60000._CUSTOM_REAL ! maxval(zstore_unique) x_target = x y_target = y ! get approximate topography elevation at target coordinates from free surface call get_topo_elevation_free_closest(x_target,y_target,elevation,distmin, & - nspec,nglob_dummy,ibool,xstore_dummy,ystore_dummy,zstore_dummy, & + nspec,nglob_unique,ibool,xstore_unique,ystore_unique,zstore_unique, & num_free_surface_faces,free_surface_ispec,free_surface_ijk) ! depth in Z-direction diff --git a/src/generate_databases/save_arrays_solver.F90 b/src/generate_databases/save_arrays_solver.F90 index 397954085..b92654e9d 100644 --- a/src/generate_databases/save_arrays_solver.F90 +++ b/src/generate_databases/save_arrays_solver.F90 @@ -60,7 +60,6 @@ subroutine save_arrays_solver_ext_mesh(nspec,ibool) nfaces_surface use create_regions_mesh_ext_par - use create_regions_mesh_ext_par, only: nglob => nglob_dummy implicit none @@ -71,10 +70,14 @@ subroutine save_arrays_solver_ext_mesh(nspec,ibool) ! local parameters integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh_dummy integer :: max_nibool_interfaces_ext_mesh + integer :: nglob integer :: ier,i,itest character(len=MAX_STRING_LEN) :: filename + ! number of unique global nodes + nglob = nglob_unique + ! database file name filename = prname(1:len_trim(prname))//'external_mesh.bin' @@ -95,9 +98,9 @@ subroutine save_arrays_solver_ext_mesh(nspec,ibool) write(IOUT) ibool - write(IOUT) xstore_dummy - write(IOUT) ystore_dummy - write(IOUT) zstore_dummy + write(IOUT) xstore_unique + write(IOUT) ystore_unique + write(IOUT) zstore_unique write(IOUT) irregular_element_number write(IOUT) xix_regular @@ -438,7 +441,6 @@ subroutine save_arrays_solver_files(nspec,ibool) use generate_databases_par, only: nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,num_interfaces_ext_mesh use create_regions_mesh_ext_par - use create_regions_mesh_ext_par, only: nglob => nglob_dummy implicit none @@ -470,19 +472,19 @@ subroutine save_arrays_solver_files(nspec,ibool) !--- x coordinate open(unit=IOUT,file=prname(1:len_trim(prname))//'x.bin',status='unknown',form='unformatted',iostat=ier) if (ier /= 0) stop 'error opening file x.bin' - write(IOUT) xstore_dummy + write(IOUT) xstore_unique close(IOUT) !--- y coordinate open(unit=IOUT,file=prname(1:len_trim(prname))//'y.bin',status='unknown',form='unformatted',iostat=ier) if (ier /= 0) stop 'error opening file y.bin' - write(IOUT) ystore_dummy + write(IOUT) ystore_unique close(IOUT) !--- z coordinate open(unit=IOUT,file=prname(1:len_trim(prname))//'z.bin',status='unknown',form='unformatted',iostat=ier) if (ier /= 0) stop 'error opening file z.bin' - write(IOUT) zstore_dummy + write(IOUT) zstore_unique close(IOUT) ! ibool @@ -511,8 +513,8 @@ subroutine save_arrays_solver_files(nspec,ibool) ! vp values - VTK file output filename = prname(1:len_trim(prname))//'vp' - call write_VTK_data_gll_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + call write_VTK_data_gll_cr(nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,ibool, & v_tmp,filename) @@ -532,8 +534,8 @@ subroutine save_arrays_solver_files(nspec,ibool) ! vs values - VTK file output filename = prname(1:len_trim(prname))//'vs' - call write_VTK_data_gll_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + call write_VTK_data_gll_cr(nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,ibool, & v_tmp,filename) ! outputs density model for check @@ -553,8 +555,8 @@ subroutine save_arrays_solver_files(nspec,ibool) ! shear attenuation - VTK file output filename = prname(1:len_trim(prname))//'qmu' - call write_VTK_data_gll_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + call write_VTK_data_gll_cr(nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,ibool, & qmu_attenuation_store,filename) ! bulk attenuation Qkappa @@ -565,8 +567,8 @@ subroutine save_arrays_solver_files(nspec,ibool) ! bulk attenuation - VTK file output filename = prname(1:len_trim(prname))//'qkappa' - call write_VTK_data_gll_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + call write_VTK_data_gll_cr(nspec,nglob_unique, & + xstore_unique,ystore_unique,zstore_unique,ibool, & qkappa_attenuation_store,filename) ! frees temporary array @@ -599,8 +601,8 @@ subroutine save_arrays_solver_files(nspec,ibool) enddo enddo filename = prname(1:len_trim(prname))//'free_surface' - call write_VTK_data_points(nglob, & - xstore_dummy,ystore_dummy,zstore_dummy, & + call write_VTK_data_points(nglob_unique, & + xstore_unique,ystore_unique,zstore_unique, & iglob_tmp,NGLLSQUARE*num_free_surface_faces, & filename) @@ -626,7 +628,7 @@ subroutine save_arrays_solver_files(nspec,ibool) enddo enddo filename = prname(1:len_trim(prname))//'coupling_acoustic_elastic' - call write_VTK_data_points(nglob,xstore_dummy,ystore_dummy,zstore_dummy, & + call write_VTK_data_points(nglob_unique,xstore_unique,ystore_unique,zstore_unique, & iglob_tmp,num_points,filename) deallocate(iglob_tmp) endif !if (ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION ) @@ -650,7 +652,7 @@ subroutine save_arrays_solver_files(nspec,ibool) enddo enddo filename = prname(1:len_trim(prname))//'coupling_acoustic_poroelastic' - call write_VTK_data_points(nglob,xstore_dummy,ystore_dummy,zstore_dummy, & + call write_VTK_data_points(nglob_unique,xstore_unique,ystore_unique,zstore_unique, & iglob_tmp,num_points,filename) deallocate(iglob_tmp) endif !if (ACOUSTIC_SIMULATION .and. POROELASTIC_SIMULATION ) @@ -674,7 +676,7 @@ subroutine save_arrays_solver_files(nspec,ibool) enddo enddo filename = prname(1:len_trim(prname))//'coupling_elastic_poroelastic' - call write_VTK_data_points(nglob,xstore_dummy,ystore_dummy,zstore_dummy, & + call write_VTK_data_points(nglob_unique,xstore_unique,ystore_unique,zstore_unique, & iglob_tmp,num_points,filename) deallocate(iglob_tmp) endif !if (ACOUSTIC_SIMULATION .and. POROELASTIC_SIMULATION @@ -699,7 +701,7 @@ subroutine save_arrays_solver_files(nspec,ibool) endif enddo filename = prname(1:len_trim(prname))//'acoustic_elastic_poroelastic_flag' - call write_VTK_data_elem_i(nspec,nglob,xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + call write_VTK_data_elem_i(nspec,nglob_unique,xstore_unique,ystore_unique,zstore_unique,ibool, & v_tmp_i,filename) deallocate(v_tmp_i) endif @@ -721,7 +723,7 @@ subroutine save_arrays_solver_files(nspec,ibool) enddo filename = prname(1:len_trim(prname))//'MPI_points' - call write_VTK_data_points(nglob,xstore_dummy,ystore_dummy,zstore_dummy, & + call write_VTK_data_points(nglob_unique,xstore_unique,ystore_unique,zstore_unique, & iglob_tmp,num_points,filename) deallocate(iglob_tmp) endif ! NPROC > 1 @@ -798,7 +800,7 @@ subroutine save_arrays_solver_injection_boundary(nspec,ibool) ny = abs_boundary_normal(2,igll,iface) nz = abs_boundary_normal(3,igll,iface) - write(IOUT,'(6f25.10)') xstore_dummy(iglob), ystore_dummy(iglob), zstore_dummy(iglob), nx, ny, nz + write(IOUT,'(6f25.10)') xstore_unique(iglob), ystore_unique(iglob), zstore_unique(iglob), nx, ny, nz enddo endif diff --git a/src/generate_databases/save_arrays_solver_adios.F90 b/src/generate_databases/save_arrays_solver_adios.F90 index ff44e5781..e592b0117 100644 --- a/src/generate_databases/save_arrays_solver_adios.F90 +++ b/src/generate_databases/save_arrays_solver_adios.F90 @@ -68,7 +68,6 @@ subroutine save_arrays_solver_ext_mesh_adios(nspec,ibool) nfaces_surface use create_regions_mesh_ext_par - use create_regions_mesh_ext_par, only: nglob => nglob_dummy use adios_helpers_mod use adios_manager_mod, only: comm_adios @@ -84,6 +83,7 @@ subroutine save_arrays_solver_ext_mesh_adios(nspec,ibool) integer :: max_nibool_interfaces_ext_mesh integer :: ier,i + integer :: nglob !--- Local parameters for ADIOS --- character(len=MAX_STRING_LEN) :: output_name @@ -116,6 +116,9 @@ subroutine save_arrays_solver_ext_mesh_adios(nspec,ibool) integer, dimension(num_vars) :: max_global_values integer :: comm + ! number of unique global nodes + nglob = nglob_unique + ! initializes output_name = LOCAL_PATH(1:len_trim(LOCAL_PATH)) // "/external_mesh.bp" @@ -268,9 +271,9 @@ subroutine save_arrays_solver_ext_mesh_adios(nspec,ibool) call define_adios_global_array1D(group, groupsize, local_dim, '', STRINGIFY_VAR(ibool)) local_dim = nglob_wmax - call define_adios_global_array1D(group, groupsize, local_dim, '', "x_global", xstore_dummy) - call define_adios_global_array1D(group, groupsize, local_dim, '', "y_global", ystore_dummy) - call define_adios_global_array1D(group, groupsize, local_dim, '', "z_global", zstore_dummy) + call define_adios_global_array1D(group, groupsize, local_dim, '', "x_global", xstore_unique) + call define_adios_global_array1D(group, groupsize, local_dim, '', "y_global", ystore_unique) + call define_adios_global_array1D(group, groupsize, local_dim, '', "z_global", zstore_unique) local_dim = nspec_wmax call define_adios_global_array1D(group, groupsize, local_dim, '', STRINGIFY_VAR(irregular_element_number)) @@ -628,9 +631,9 @@ subroutine save_arrays_solver_ext_mesh_adios(nspec,ibool) call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, STRINGIFY_VAR(ibool)) local_dim = nglob_wmax - call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "x_global", xstore_dummy) - call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "y_global", ystore_dummy) - call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "z_global", zstore_dummy) + call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "x_global", xstore_unique) + call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "y_global", ystore_unique) + call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "z_global", zstore_unique) local_dim = nspec_wmax call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, STRINGIFY_VAR(irregular_element_number)) @@ -991,7 +994,7 @@ subroutine save_arrays_solver_ext_mesh_adios(nspec,ibool) !if (num_interfaces_ext_mesh >= 1) then ! filename = prname(1:len_trim(prname))//'MPI_1_points' ! call write_VTK_data_points(nglob, & - ! xstore_dummy,ystore_dummy,zstore_dummy, & + ! xstore_unique,ystore_unique,zstore_unique, & ! ibool_interfaces_ext_mesh_dummy(1:nibool_interfaces_ext_mesh(1),1), & ! nibool_interfaces_ext_mesh(1), & ! filename) @@ -1110,9 +1113,9 @@ subroutine save_arrays_solver_files_adios(nspec, nglob, ibool, nspec_wmax, nglob call define_adios_scalar(group, groupsize, '', STRINGIFY_VAR(nglob)) local_dim = nglob_wmax - call define_adios_global_array1D(group, groupsize, local_dim, '', "x_global", xstore_dummy) - call define_adios_global_array1D(group, groupsize, local_dim, '', "y_global", ystore_dummy) - call define_adios_global_array1D(group, groupsize, local_dim, '', "z_global", zstore_dummy) + call define_adios_global_array1D(group, groupsize, local_dim, '', "x_global", xstore_unique) + call define_adios_global_array1D(group, groupsize, local_dim, '', "y_global", ystore_unique) + call define_adios_global_array1D(group, groupsize, local_dim, '', "z_global", zstore_unique) local_dim = NGLLX * NGLLY * NGLLZ * nspec_wmax call define_adios_global_array1D(group, groupsize, local_dim, '', STRINGIFY_VAR(xstore)) @@ -1145,9 +1148,9 @@ subroutine save_arrays_solver_files_adios(nspec, nglob, ibool, nspec_wmax, nglob call check_adios_err(myrank,ier) local_dim = nglob_wmax - call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "x_global", xstore_dummy) - call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "y_global", ystore_dummy) - call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "z_global", zstore_dummy) + call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "x_global", xstore_unique) + call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "y_global", ystore_unique) + call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, "z_global", zstore_unique) local_dim = NGLLX * NGLLY * NGLLZ * nspec_wmax call write_adios_global_1d_array(handle, myrank, sizeprocs, local_dim, STRINGIFY_VAR(xstore)) diff --git a/src/generate_databases/setup_color_perm.f90 b/src/generate_databases/setup_color_perm.f90 index 89045ab99..a6c5e90be 100644 --- a/src/generate_databases/setup_color_perm.f90 +++ b/src/generate_databases/setup_color_perm.f90 @@ -245,7 +245,7 @@ subroutine setup_color(myrank,nspec,nglob,ibool,perm,ispec_is_d,idomain, & if (SAVE_MESH_FILES .and. DEBUG) then filename = prname(1:len_trim(prname))//'color_'//str_domain(idomain) call write_VTK_data_elem_i(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore_unique,ystore_unique,zstore_unique,ibool, & color,filename) endif @@ -403,7 +403,7 @@ subroutine setup_color(myrank,nspec,nglob,ibool,perm,ispec_is_d,idomain, & ! debug: outputs permutation array as vtk file if (DEBUG) then filename = prname(1:len_trim(prname))//'perm_'//str_domain(idomain) - call write_VTK_data_elem_i(nspec,nglob,xstore_dummy,ystore_dummy,zstore_dummy,ibool,perm,filename) + call write_VTK_data_elem_i(nspec,nglob,xstore_unique,ystore_unique,zstore_unique,ibool,perm,filename) endif deallocate(num_of_elems_in_this_color) @@ -573,7 +573,7 @@ subroutine setup_permutation(myrank,nspec,nglob,ibool,ANISOTROPY,perm,SAVE_MESH_ if (SAVE_MESH_FILES .and. DEBUG) then filename = prname(1:len_trim(prname))//'perm_global' call write_VTK_data_elem_i(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore_unique,ystore_unique,zstore_unique,ibool, & temp_perm_global,filename) endif diff --git a/src/generate_databases/setup_mesh.f90 b/src/generate_databases/setup_mesh.f90 index d12973788..ff52b2104 100644 --- a/src/generate_databases/setup_mesh.f90 +++ b/src/generate_databases/setup_mesh.f90 @@ -41,7 +41,7 @@ subroutine setup_mesh ! compute maximum number of points npointot = NSPEC_AB * NGLLX * NGLLY * NGLLZ -! use dynamic allocation to allocate memory for arrays + ! use dynamic allocation to allocate memory for arrays allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 605') if (ier /= 0) stop 'error allocating array ibool' @@ -63,16 +63,16 @@ subroutine setup_mesh max_memory_size = max_memory_size_request -! make sure everybody is synchronized + ! make sure everybody is synchronized call synchronize_all() -! main working routine to create all the regions of the mesh + ! main working routine to create all the regions of the mesh if (myrank == 0) then write(IMAIN,*) 'create regions:' endif call create_regions_mesh() -! print min and max of topography included + ! print min and max of topography included min_elevation = HUGEVAL max_elevation = -HUGEVAL do iface = 1,nspec2D_top_ext @@ -87,7 +87,7 @@ subroutine setup_mesh enddo enddo -! compute the maximum of the maxima for all the slices using an MPI reduction + ! compute the maximum of the maxima for all the slices using an MPI reduction call min_all_dp(min_elevation,min_elevation_all) call max_all_dp(max_elevation,max_elevation_all) @@ -106,13 +106,22 @@ subroutine setup_mesh call flush_IMAIN() endif -! clean-up + ! clean-up deallocate(xstore,ystore,zstore) deallocate(ibool) deallocate(ispec_is_surface_external_mesh) deallocate(iglob_is_surface_external_mesh) -! make sure everybody is synchronized + ! make sure everybody is synchronized + call synchronize_all() + + ! user output + if (myrank == 0) then + write(IMAIN,*) + write(IMAIN,*) 'done mesh setup' + write(IMAIN,*) + call flush_IMAIN() + endif call synchronize_all() end subroutine setup_mesh diff --git a/src/inverse_problem_for_model/input_output/Teleseismic_IO.f90 b/src/inverse_problem_for_model/input_output/Teleseismic_IO.f90 index e6953cd7d..7fcb69443 100644 --- a/src/inverse_problem_for_model/input_output/Teleseismic_IO.f90 +++ b/src/inverse_problem_for_model/input_output/Teleseismic_IO.f90 @@ -468,9 +468,9 @@ subroutine setup_teleseismic_stations(acqui_simu, myrank) ! get mesh properties (mandatory before calling locate_point_in_mesh) call usual_hex_nodes(NGNOD,iaddx,iaddy,iaddz) call check_mesh_distances(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & - x_min_glob,x_max_glob,y_min_glob,y_max_glob,z_min_glob,z_max_glob, & - elemsize_min_glob,elemsize_max_glob, & - distance_min_glob,distance_max_glob) + x_min_glob,x_max_glob,y_min_glob,y_max_glob,z_min_glob,z_max_glob, & + elemsize_min_glob,elemsize_max_glob, & + distance_min_glob,distance_max_glob) if (DEBUG_MODE) then write (IIDD, *) write (IIDD, *) ' LOCATE TELESEISMIC STATIONS ---------------------' diff --git a/src/inverse_problem_for_model/input_output/input_output_mod.f90 b/src/inverse_problem_for_model/input_output/input_output_mod.f90 index 11ace429e..ffac89ca9 100644 --- a/src/inverse_problem_for_model/input_output/input_output_mod.f90 +++ b/src/inverse_problem_for_model/input_output/input_output_mod.f90 @@ -238,12 +238,12 @@ subroutine init_input_output_mod(inversion_param, acqui_simu, myrank) !! get dimension of mesh call check_mesh_distances(myrank,NSPEC_AB,NGLOB_AB, & - ibool,xstore,ystore,zstore, & - inversion_param%xmin,inversion_param%xmax, & - inversion_param%ymin,inversion_param%ymax, & - inversion_param%zmin,inversion_param%zmax, & - elemsize_min_glob,elemsize_max_glob, & - distance_min_glob,distance_max_glob) + ibool,xstore,ystore,zstore, & + inversion_param%xmin,inversion_param%xmax, & + inversion_param%ymin,inversion_param%ymax, & + inversion_param%zmin,inversion_param%zmax, & + elemsize_min_glob,elemsize_max_glob, & + distance_min_glob,distance_max_glob) !!------------------------------------------------------------------------------------------------------------------------- diff --git a/src/inverse_problem_for_model/projection_on_FD_grid/projection_on_FD_grid_mod.f90 b/src/inverse_problem_for_model/projection_on_FD_grid/projection_on_FD_grid_mod.f90 index cc681b161..f364168fc 100644 --- a/src/inverse_problem_for_model/projection_on_FD_grid/projection_on_FD_grid_mod.f90 +++ b/src/inverse_problem_for_model/projection_on_FD_grid/projection_on_FD_grid_mod.f90 @@ -193,11 +193,11 @@ subroutine compute_interpolation_coeff_FD_SEM(projection_fd, myrank) ! get mesh properties (mandatory before calling locate_point_in_mesh_simple) call usual_hex_nodes(NGNOD,iaddx,iaddy,iaddz) call check_mesh_distances(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & - x_min_glob,x_max_glob,y_min_glob,y_max_glob,z_min_glob,z_max_glob, & - elemsize_min_glob,elemsize_max_glob, & - distance_min_glob,distance_max_glob) + x_min_glob,x_max_glob,y_min_glob,y_max_glob,z_min_glob,z_max_glob, & + elemsize_min_glob,elemsize_max_glob, & + distance_min_glob,distance_max_glob) - nb_fd_point_loc=0 + nb_fd_point_loc = 0 allocate(point_already_found(nx_fd_proj, ny_fd_proj, nz_fd_proj),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 144') @@ -208,7 +208,7 @@ subroutine compute_interpolation_coeff_FD_SEM(projection_fd, myrank) allocate(gamma_in_fd(nx_fd_proj, ny_fd_proj, nz_fd_proj),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 147') - point_already_found(:,:,:)=0 + point_already_found(:,:,:) = 0 !! loop over elements do ispec = 1, NSPEC_AB diff --git a/src/inverse_problem_for_model/regularization/regularization_SEM_mod.f90 b/src/inverse_problem_for_model/regularization/regularization_SEM_mod.f90 index 7820cf698..a0ea25a9f 100644 --- a/src/inverse_problem_for_model/regularization/regularization_SEM_mod.f90 +++ b/src/inverse_problem_for_model/regularization/regularization_SEM_mod.f90 @@ -2373,29 +2373,27 @@ subroutine compute_mean_values_on_edge(field) enddo enddo - field_wksp(:,:) = field_wksp(:,:) / valence(:,:) valence(:,:) = 1. !! 2/ compute average values of field_to_derivate at the edge of MPI slices ------------------- call assemble_MPI_vector_blocking(NPROC,NGLOB_AB,valence, & - num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & - nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & - my_neighbors_ext_mesh_sp) + num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & + nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & + my_neighbors_ext_mesh_sp) call assemble_MPI_vector_blocking(NPROC,NGLOB_AB,field_wksp, & - num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & - nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & - my_neighbors_ext_mesh_sp) + num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & + nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & + my_neighbors_ext_mesh_sp) field_wksp(:,:) = field_wksp(:,:) / valence(:,:) - - do ispec=1, NSPEC_AB - do k=1,NGLLZ - do j=1,NGLLY - do i=1,NGLLX + do ispec = 1, NSPEC_AB + do k = 1,NGLLZ + do j = 1,NGLLY + do i = 1,NGLLX iglob = ibool(i,j,k,ispec) field(1,i,j,k,ispec) = field_wksp(1,iglob) @@ -4065,25 +4063,24 @@ subroutine compute_laplacian_FD(Dfb, field_input, regularization_fd) if (ier /= 0) call exit_MPI_without_rank('error allocating array 200') allocate(field_to_send(NGLOB_AB), field_overlap(indx_recv(NPROC)), result_df(nline),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 201') - valence(:,:)=1. + valence(:,:) = 1. field_to_derivate(:,:) = field_input(:,:) !! 1/ compute average values of field_to_derivate at the edge of MPI slices ------------------- call assemble_MPI_vector_blocking(NPROC,NGLOB_AB,valence, & - num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & - nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & - my_neighbors_ext_mesh_sp) + num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & + nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & + my_neighbors_ext_mesh_sp) call assemble_MPI_vector_blocking(NPROC,NGLOB_AB,field_to_derivate, & - num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & - nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & - my_neighbors_ext_mesh_sp) + num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & + nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & + my_neighbors_ext_mesh_sp) !! divide by valence field_to_derivate(:,:) = field_to_derivate(:,:) / valence(:,:) - - do idim=1, NDIM + do idim = 1, NDIM !! 2/ communicate overlap of MPI slices field_to_send(:)=field_to_derivate(idim,:) call send_recv_blocking(field_to_send, field_overlap) @@ -4149,7 +4146,7 @@ subroutine compute_gradient_laplacian_FD(nGrad, Laplac, field_input, regularizat if (ier /= 0) call exit_MPI_without_rank('error allocating array 203') allocate(field_to_send(NGLOB_AB), field_overlap(indx_recv(NPROC)), result_df(nline),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 204') - valence(:,:)=1. + valence(:,:) = 1. !! need to duplicate in order to use the already build subroutine from sepcfem package : assembel_MPI.. field_to_derivate(1,:) = field_input(:) @@ -4158,14 +4155,14 @@ subroutine compute_gradient_laplacian_FD(nGrad, Laplac, field_input, regularizat !! 1/ compute average values of field_to_derivate at the edge of MPI slices ------------------- call assemble_MPI_vector_blocking(NPROC,NGLOB_AB,valence, & - num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & - nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & - my_neighbors_ext_mesh_sp) + num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & + nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & + my_neighbors_ext_mesh_sp) call assemble_MPI_vector_blocking(NPROC,NGLOB_AB,field_to_derivate, & - num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & - nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & - my_neighbors_ext_mesh_sp) + num_interfaces_ext_mesh_sp,max_nibool_interfaces_ext_mesh_sp, & + nibool_interfaces_ext_mesh_sp,ibool_interfaces_ext_mesh_sp, & + my_neighbors_ext_mesh_sp) !! divide by valence field_to_derivate(1,:) = field_to_derivate(1,:) / valence(1,:) diff --git a/src/shared/get_element_face.f90 b/src/shared/get_element_face.f90 index 4f62af00a..080cc2aac 100644 --- a/src/shared/get_element_face.f90 +++ b/src/shared/get_element_face.f90 @@ -27,7 +27,7 @@ subroutine get_element_face_id(ispec,xcoord,ycoord,zcoord, & ibool,nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy, & + xstore_unique,ystore_unique,zstore_unique, & iface_id ) ! returns iface_id of face in reference element, determined by corner locations xcoord/ycoord/zcoord; @@ -46,7 +46,7 @@ subroutine get_element_face_id(ispec,xcoord,ycoord,zcoord, & integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool ! global point locations - real(kind=CUSTOM_REAL),dimension(nglob),intent(in) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL),dimension(nglob),intent(in) :: xstore_unique,ystore_unique,zstore_unique ! local parameters real(kind=CUSTOM_REAL) :: xcoord_face,ycoord_face,zcoord_face @@ -109,9 +109,9 @@ subroutine get_element_face_id(ispec,xcoord,ycoord,zcoord, & ! coordinates iglob = ibool(i,j,k,ispec) - xcoord_face = xstore_dummy(iglob) - ycoord_face = ystore_dummy(iglob) - zcoord_face = zstore_dummy(iglob) + xcoord_face = xstore_unique(iglob) + ycoord_face = ystore_unique(iglob) + zcoord_face = zstore_unique(iglob) ! face midpoint coordinates midpoint_faces(1,ifa) = midpoint_faces(1,ifa) + xcoord_face * avg_factor @@ -143,7 +143,7 @@ subroutine get_element_face_id(ispec,xcoord,ycoord,zcoord, & j = iface_all_corner_ijk(2,icorner,iloc(1)) k = iface_all_corner_ijk(3,icorner,iloc(1)) iglob = ibool(i,j,k,ispec) - print *,'error corner:',icorner,'xyz:',xstore_dummy(iglob),ystore_dummy(iglob),zstore_dummy(iglob) + print *,'error corner:',icorner,'xyz:',xstore_unique(iglob),ystore_unique(iglob),zstore_unique(iglob) enddo ! target do icorner = 1,NGNOD2D_FOUR_CORNERS @@ -279,7 +279,7 @@ end subroutine get_element_face_gll_indices subroutine get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & ibool,nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy, & + xstore_unique,ystore_unique,zstore_unique, & normal) ! only changes direction of normal to point outwards of element @@ -305,7 +305,7 @@ subroutine get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool ! global point locations - real(kind=CUSTOM_REAL),dimension(nglob),intent(in) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL),dimension(nglob),intent(in) :: xstore_unique,ystore_unique,zstore_unique ! face normal real(kind=CUSTOM_REAL),dimension(NDIM),intent(inout) :: normal @@ -349,9 +349,9 @@ subroutine get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, & end select ! vector from corner 1 to this opposite one - v_tmp(1) = xstore_dummy(iglob) - xcoord(1) - v_tmp(2) = ystore_dummy(iglob) - ycoord(1) - v_tmp(3) = zstore_dummy(iglob) - zcoord(1) + v_tmp(1) = xstore_unique(iglob) - xcoord(1) + v_tmp(2) = ystore_unique(iglob) - ycoord(1) + v_tmp(3) = zstore_unique(iglob) - zcoord(1) ! scalar product (dot product) tmp = v_tmp(1)*face_n(1) + v_tmp(2)*face_n(2) + v_tmp(3)*face_n(3) @@ -386,7 +386,7 @@ end subroutine get_element_face_normal subroutine get_element_face_normal_idirect(ispec,iface,xcoord,ycoord,zcoord, & ibool,nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy, & + xstore_unique,ystore_unique,zstore_unique, & normal,idirect) ! returns direction of normal: @@ -414,7 +414,7 @@ subroutine get_element_face_normal_idirect(ispec,iface,xcoord,ycoord,zcoord, & integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool ! global point locations - real(kind=CUSTOM_REAL),dimension(nglob),intent(in) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL),dimension(nglob),intent(in) :: xstore_unique,ystore_unique,zstore_unique ! face normal real(kind=CUSTOM_REAL),dimension(NDIM),intent(inout) :: normal @@ -461,9 +461,9 @@ subroutine get_element_face_normal_idirect(ispec,iface,xcoord,ycoord,zcoord, & end select ! vector from corner 1 to this opposite one - v_tmp(1) = xstore_dummy(iglob) - xcoord(1) - v_tmp(2) = ystore_dummy(iglob) - ycoord(1) - v_tmp(3) = zstore_dummy(iglob) - zcoord(1) + v_tmp(1) = xstore_unique(iglob) - xcoord(1) + v_tmp(2) = ystore_unique(iglob) - ycoord(1) + v_tmp(3) = zstore_unique(iglob) - zcoord(1) ! scalar product (dot product) tmp = v_tmp(1)*face_n(1) + v_tmp(2)*face_n(2) + v_tmp(3)*face_n(3) @@ -497,7 +497,7 @@ end subroutine get_element_face_normal_idirect ! subroutine get_element_corners(ispec,iface_ref,xcoord,ycoord,zcoord,iglob_corners_ref, & - ibool,nspec,nglob,xstore_dummy,ystore_dummy,zstore_dummy, & + ibool,nspec,nglob,xstore_unique,ystore_unique,zstore_unique, & iface_all_corner_ijk) use constants, only: NGNOD2D_FOUR_CORNERS,CUSTOM_REAL,NGLLX,NGLLY,NGLLZ @@ -514,7 +514,7 @@ subroutine get_element_corners(ispec,iface_ref,xcoord,ycoord,zcoord,iglob_corner ! index array integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool ! global point locations - real(kind=CUSTOM_REAL),dimension(nglob),intent(in) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL),dimension(nglob),intent(in) :: xstore_unique,ystore_unique,zstore_unique ! assumes NGNOD2D_FOUR_CORNERS == 4 !! DK DK Oct 2012: in principle we should use NGNOD2D instead of NGNOD2D_FOUR_CORNERS when @@ -540,9 +540,9 @@ subroutine get_element_corners(ispec,iface_ref,xcoord,ycoord,zcoord,iglob_corner iglob_corners_ref(icorner) = ibool(i,j,k,ispec) ! reference corner coordinates - xcoord(icorner) = xstore_dummy(iglob_corners_ref(icorner)) - ycoord(icorner) = ystore_dummy(iglob_corners_ref(icorner)) - zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner)) + xcoord(icorner) = xstore_unique(iglob_corners_ref(icorner)) + ycoord(icorner) = ystore_unique(iglob_corners_ref(icorner)) + zcoord(icorner) = zstore_unique(iglob_corners_ref(icorner)) enddo end subroutine get_element_corners diff --git a/src/shared/get_jacobian_boundaries.f90 b/src/shared/get_jacobian_boundaries.f90 index f65d027dc..846d60b59 100644 --- a/src/shared/get_jacobian_boundaries.f90 +++ b/src/shared/get_jacobian_boundaries.f90 @@ -26,189 +26,189 @@ !===================================================================== subroutine get_jacobian_boundary_face(nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, & + xstore,ystore,zstore,ibool,nglob, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & ispec,iface,jacobian2Dw_face,normal_face,NGLLA,NGLLB,NGNOD2D) ! returns jacobian2Dw_face and normal_face (pointing outwards of element) - use constants + use constants, only: CUSTOM_REAL,NDIM,NGLLX,NGLLY,NGLLZ,MIDX,MIDY,MIDZ,NDIM2D,myrank implicit none - integer :: nspec,nglob + integer,intent(in) :: nspec,nglob -! arrays with the mesh - integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool - real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + ! arrays with the mesh + integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool + real(kind=CUSTOM_REAL), dimension(nglob),intent(in) :: xstore,ystore,zstore -! face information - integer :: iface,ispec,NGLLA,NGLLB - real(kind=CUSTOM_REAL), dimension(NGLLA,NGLLB) :: jacobian2Dw_face - real(kind=CUSTOM_REAL), dimension(NDIM,NGLLA,NGLLB) :: normal_face + ! face information + integer,intent(in) :: iface,ispec,NGLLA,NGLLB + real(kind=CUSTOM_REAL), dimension(NGLLA,NGLLB),intent(inout) :: jacobian2Dw_face + real(kind=CUSTOM_REAL), dimension(NDIM,NGLLA,NGLLB),intent(inout) :: normal_face - integer :: NGNOD2D - double precision dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ) - double precision dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ) - double precision dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY) - double precision dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY) + integer,intent(in) :: NGNOD2D + double precision,intent(in) :: dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ) + double precision,intent(in) :: dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ) + double precision,intent(in) :: dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY) + double precision,intent(in) :: dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY) - double precision, dimension(NGLLX,NGLLY) :: wgllwgll_xy - double precision, dimension(NGLLX,NGLLZ) :: wgllwgll_xz - double precision, dimension(NGLLY,NGLLZ) :: wgllwgll_yz + double precision, dimension(NGLLX,NGLLY),intent(in) :: wgllwgll_xy + double precision, dimension(NGLLX,NGLLZ),intent(in) :: wgllwgll_xz + double precision, dimension(NGLLY,NGLLZ),intent(in) :: wgllwgll_yz ! local parameters -! face corners + ! face corners double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D) -! check that the parameter file is correct + ! check that the parameter file is correct if (NGNOD2D /= 4 .and. NGNOD2D /= 9) & call exit_MPI(myrank,'surface elements should have 4 or 9 control nodes') select case (iface) ! on reference face: xmin case (1) - xelm(1)=xstore_dummy( ibool(1,1,1,ispec) ) - yelm(1)=ystore_dummy( ibool(1,1,1,ispec) ) - zelm(1)=zstore_dummy( ibool(1,1,1,ispec) ) - xelm(2)=xstore_dummy( ibool(1,NGLLY,1,ispec) ) - yelm(2)=ystore_dummy( ibool(1,NGLLY,1,ispec) ) - zelm(2)=zstore_dummy( ibool(1,NGLLY,1,ispec) ) - xelm(3)=xstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) ) - yelm(3)=ystore_dummy( ibool(1,NGLLY,NGLLZ,ispec) ) - zelm(3)=zstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) ) - xelm(4)=xstore_dummy( ibool(1,1,NGLLZ,ispec) ) - yelm(4)=ystore_dummy( ibool(1,1,NGLLZ,ispec) ) - zelm(4)=zstore_dummy( ibool(1,1,NGLLZ,ispec) ) + xelm(1)=xstore( ibool(1,1,1,ispec) ) + yelm(1)=ystore( ibool(1,1,1,ispec) ) + zelm(1)=zstore( ibool(1,1,1,ispec) ) + xelm(2)=xstore( ibool(1,NGLLY,1,ispec) ) + yelm(2)=ystore( ibool(1,NGLLY,1,ispec) ) + zelm(2)=zstore( ibool(1,NGLLY,1,ispec) ) + xelm(3)=xstore( ibool(1,NGLLY,NGLLZ,ispec) ) + yelm(3)=ystore( ibool(1,NGLLY,NGLLZ,ispec) ) + zelm(3)=zstore( ibool(1,NGLLY,NGLLZ,ispec) ) + xelm(4)=xstore( ibool(1,1,NGLLZ,ispec) ) + yelm(4)=ystore( ibool(1,1,NGLLZ,ispec) ) + zelm(4)=zstore( ibool(1,1,NGLLZ,ispec) ) if (NGNOD2D == 9) then - xelm(5)=xstore_dummy( ibool(1,MIDY,1,ispec) ) - yelm(5)=ystore_dummy( ibool(1,MIDY,1,ispec) ) - zelm(5)=zstore_dummy( ibool(1,MIDY,1,ispec) ) - xelm(6)=xstore_dummy( ibool(1,NGLLY,MIDZ,ispec) ) - yelm(6)=ystore_dummy( ibool(1,NGLLY,MIDZ,ispec) ) - zelm(6)=zstore_dummy( ibool(1,NGLLY,MIDZ,ispec) ) - xelm(7)=xstore_dummy( ibool(1,MIDY,NGLLZ,ispec) ) - yelm(7)=ystore_dummy( ibool(1,MIDY,NGLLZ,ispec) ) - zelm(7)=zstore_dummy( ibool(1,MIDY,NGLLZ,ispec) ) - xelm(8)=xstore_dummy( ibool(1,1,MIDZ,ispec) ) - yelm(8)=ystore_dummy( ibool(1,1,MIDZ,ispec) ) - zelm(8)=zstore_dummy( ibool(1,1,MIDZ,ispec) ) - xelm(9)=xstore_dummy( ibool(1,MIDY,MIDZ,ispec) ) - yelm(9)=ystore_dummy( ibool(1,MIDY,MIDZ,ispec) ) - zelm(9)=zstore_dummy( ibool(1,MIDY,MIDZ,ispec) ) + xelm(5)=xstore( ibool(1,MIDY,1,ispec) ) + yelm(5)=ystore( ibool(1,MIDY,1,ispec) ) + zelm(5)=zstore( ibool(1,MIDY,1,ispec) ) + xelm(6)=xstore( ibool(1,NGLLY,MIDZ,ispec) ) + yelm(6)=ystore( ibool(1,NGLLY,MIDZ,ispec) ) + zelm(6)=zstore( ibool(1,NGLLY,MIDZ,ispec) ) + xelm(7)=xstore( ibool(1,MIDY,NGLLZ,ispec) ) + yelm(7)=ystore( ibool(1,MIDY,NGLLZ,ispec) ) + zelm(7)=zstore( ibool(1,MIDY,NGLLZ,ispec) ) + xelm(8)=xstore( ibool(1,1,MIDZ,ispec) ) + yelm(8)=ystore( ibool(1,1,MIDZ,ispec) ) + zelm(8)=zstore( ibool(1,1,MIDZ,ispec) ) + xelm(9)=xstore( ibool(1,MIDY,MIDZ,ispec) ) + yelm(9)=ystore( ibool(1,MIDY,MIDZ,ispec) ) + zelm(9)=zstore( ibool(1,MIDY,MIDZ,ispec) ) endif call compute_jacobian_2D_face(xelm,yelm,zelm, & dershape2D_x,wgllwgll_yz, & jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) -! on boundary: xmax + ! on boundary: xmax case (2) - xelm(1)=xstore_dummy( ibool(NGLLX,1,1,ispec) ) - yelm(1)=ystore_dummy( ibool(NGLLX,1,1,ispec) ) - zelm(1)=zstore_dummy( ibool(NGLLX,1,1,ispec) ) - xelm(2)=xstore_dummy( ibool(NGLLX,NGLLY,1,ispec) ) - yelm(2)=ystore_dummy( ibool(NGLLX,NGLLY,1,ispec) ) - zelm(2)=zstore_dummy( ibool(NGLLX,NGLLY,1,ispec) ) - xelm(3)=xstore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) - yelm(3)=ystore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) - zelm(3)=zstore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) - xelm(4)=xstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) ) - yelm(4)=ystore_dummy( ibool(NGLLX,1,NGLLZ,ispec) ) - zelm(4)=zstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) ) + xelm(1)=xstore( ibool(NGLLX,1,1,ispec) ) + yelm(1)=ystore( ibool(NGLLX,1,1,ispec) ) + zelm(1)=zstore( ibool(NGLLX,1,1,ispec) ) + xelm(2)=xstore( ibool(NGLLX,NGLLY,1,ispec) ) + yelm(2)=ystore( ibool(NGLLX,NGLLY,1,ispec) ) + zelm(2)=zstore( ibool(NGLLX,NGLLY,1,ispec) ) + xelm(3)=xstore( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) + yelm(3)=ystore( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) + zelm(3)=zstore( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) + xelm(4)=xstore( ibool(NGLLX,1,NGLLZ,ispec) ) + yelm(4)=ystore( ibool(NGLLX,1,NGLLZ,ispec) ) + zelm(4)=zstore( ibool(NGLLX,1,NGLLZ,ispec) ) if (NGNOD2D == 9) then - xelm(5)=xstore_dummy( ibool(NGLLX,MIDY,1,ispec) ) - yelm(5)=ystore_dummy( ibool(NGLLX,MIDY,1,ispec) ) - zelm(5)=zstore_dummy( ibool(NGLLX,MIDY,1,ispec) ) - xelm(6)=xstore_dummy( ibool(NGLLX,NGLLY,MIDZ,ispec) ) - yelm(6)=ystore_dummy( ibool(NGLLX,NGLLY,MIDZ,ispec) ) - zelm(6)=zstore_dummy( ibool(NGLLX,NGLLY,MIDZ,ispec) ) - xelm(7)=xstore_dummy( ibool(NGLLX,MIDY,NGLLZ,ispec) ) - yelm(7)=ystore_dummy( ibool(NGLLX,MIDY,NGLLZ,ispec) ) - zelm(7)=zstore_dummy( ibool(NGLLX,MIDY,NGLLZ,ispec) ) - xelm(8)=xstore_dummy( ibool(NGLLX,1,MIDZ,ispec) ) - yelm(8)=ystore_dummy( ibool(NGLLX,1,MIDZ,ispec) ) - zelm(8)=zstore_dummy( ibool(NGLLX,1,MIDZ,ispec) ) - xelm(9)=xstore_dummy( ibool(NGLLX,MIDY,MIDZ,ispec) ) - yelm(9)=ystore_dummy( ibool(NGLLX,MIDY,MIDZ,ispec) ) - zelm(9)=zstore_dummy( ibool(NGLLX,MIDY,MIDZ,ispec) ) + xelm(5)=xstore( ibool(NGLLX,MIDY,1,ispec) ) + yelm(5)=ystore( ibool(NGLLX,MIDY,1,ispec) ) + zelm(5)=zstore( ibool(NGLLX,MIDY,1,ispec) ) + xelm(6)=xstore( ibool(NGLLX,NGLLY,MIDZ,ispec) ) + yelm(6)=ystore( ibool(NGLLX,NGLLY,MIDZ,ispec) ) + zelm(6)=zstore( ibool(NGLLX,NGLLY,MIDZ,ispec) ) + xelm(7)=xstore( ibool(NGLLX,MIDY,NGLLZ,ispec) ) + yelm(7)=ystore( ibool(NGLLX,MIDY,NGLLZ,ispec) ) + zelm(7)=zstore( ibool(NGLLX,MIDY,NGLLZ,ispec) ) + xelm(8)=xstore( ibool(NGLLX,1,MIDZ,ispec) ) + yelm(8)=ystore( ibool(NGLLX,1,MIDZ,ispec) ) + zelm(8)=zstore( ibool(NGLLX,1,MIDZ,ispec) ) + xelm(9)=xstore( ibool(NGLLX,MIDY,MIDZ,ispec) ) + yelm(9)=ystore( ibool(NGLLX,MIDY,MIDZ,ispec) ) + zelm(9)=zstore( ibool(NGLLX,MIDY,MIDZ,ispec) ) endif call compute_jacobian_2D_face(xelm,yelm,zelm, & dershape2D_x,wgllwgll_yz, & jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) -! on boundary: ymin + ! on boundary: ymin case (3) - xelm(1)=xstore_dummy( ibool(1,1,1,ispec) ) - yelm(1)=ystore_dummy( ibool(1,1,1,ispec) ) - zelm(1)=zstore_dummy( ibool(1,1,1,ispec) ) - xelm(2)=xstore_dummy( ibool(NGLLX,1,1,ispec) ) - yelm(2)=ystore_dummy( ibool(NGLLX,1,1,ispec) ) - zelm(2)=zstore_dummy( ibool(NGLLX,1,1,ispec) ) - xelm(3)=xstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) ) - yelm(3)=ystore_dummy( ibool(NGLLX,1,NGLLZ,ispec) ) - zelm(3)=zstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) ) - xelm(4)=xstore_dummy( ibool(1,1,NGLLZ,ispec) ) - yelm(4)=ystore_dummy( ibool(1,1,NGLLZ,ispec) ) - zelm(4)=zstore_dummy( ibool(1,1,NGLLZ,ispec) ) + xelm(1)=xstore( ibool(1,1,1,ispec) ) + yelm(1)=ystore( ibool(1,1,1,ispec) ) + zelm(1)=zstore( ibool(1,1,1,ispec) ) + xelm(2)=xstore( ibool(NGLLX,1,1,ispec) ) + yelm(2)=ystore( ibool(NGLLX,1,1,ispec) ) + zelm(2)=zstore( ibool(NGLLX,1,1,ispec) ) + xelm(3)=xstore( ibool(NGLLX,1,NGLLZ,ispec) ) + yelm(3)=ystore( ibool(NGLLX,1,NGLLZ,ispec) ) + zelm(3)=zstore( ibool(NGLLX,1,NGLLZ,ispec) ) + xelm(4)=xstore( ibool(1,1,NGLLZ,ispec) ) + yelm(4)=ystore( ibool(1,1,NGLLZ,ispec) ) + zelm(4)=zstore( ibool(1,1,NGLLZ,ispec) ) if (NGNOD2D == 9) then - xelm(5)=xstore_dummy( ibool(MIDX,1,1,ispec) ) - yelm(5)=ystore_dummy( ibool(MIDX,1,1,ispec) ) - zelm(5)=zstore_dummy( ibool(MIDX,1,1,ispec) ) - xelm(6)=xstore_dummy( ibool(NGLLX,1,MIDZ,ispec) ) - yelm(6)=ystore_dummy( ibool(NGLLX,1,MIDZ,ispec) ) - zelm(6)=zstore_dummy( ibool(NGLLX,1,MIDZ,ispec) ) - xelm(7)=xstore_dummy( ibool(MIDX,1,NGLLZ,ispec) ) - yelm(7)=ystore_dummy( ibool(MIDX,1,NGLLZ,ispec) ) - zelm(7)=zstore_dummy( ibool(MIDX,1,NGLLZ,ispec) ) - xelm(8)=xstore_dummy( ibool(1,1,MIDZ,ispec) ) - yelm(8)=ystore_dummy( ibool(1,1,MIDZ,ispec) ) - zelm(8)=zstore_dummy( ibool(1,1,MIDZ,ispec) ) - xelm(9)=xstore_dummy( ibool(MIDX,1,MIDZ,ispec) ) - yelm(9)=ystore_dummy( ibool(MIDX,1,MIDZ,ispec) ) - zelm(9)=zstore_dummy( ibool(MIDX,1,MIDZ,ispec) ) + xelm(5)=xstore( ibool(MIDX,1,1,ispec) ) + yelm(5)=ystore( ibool(MIDX,1,1,ispec) ) + zelm(5)=zstore( ibool(MIDX,1,1,ispec) ) + xelm(6)=xstore( ibool(NGLLX,1,MIDZ,ispec) ) + yelm(6)=ystore( ibool(NGLLX,1,MIDZ,ispec) ) + zelm(6)=zstore( ibool(NGLLX,1,MIDZ,ispec) ) + xelm(7)=xstore( ibool(MIDX,1,NGLLZ,ispec) ) + yelm(7)=ystore( ibool(MIDX,1,NGLLZ,ispec) ) + zelm(7)=zstore( ibool(MIDX,1,NGLLZ,ispec) ) + xelm(8)=xstore( ibool(1,1,MIDZ,ispec) ) + yelm(8)=ystore( ibool(1,1,MIDZ,ispec) ) + zelm(8)=zstore( ibool(1,1,MIDZ,ispec) ) + xelm(9)=xstore( ibool(MIDX,1,MIDZ,ispec) ) + yelm(9)=ystore( ibool(MIDX,1,MIDZ,ispec) ) + zelm(9)=zstore( ibool(MIDX,1,MIDZ,ispec) ) endif call compute_jacobian_2D_face(xelm,yelm,zelm, & dershape2D_y,wgllwgll_xz, & jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) -! on boundary: ymax + ! on boundary: ymax case (4) - xelm(1)=xstore_dummy( ibool(1,NGLLY,1,ispec) ) - yelm(1)=ystore_dummy( ibool(1,NGLLY,1,ispec) ) - zelm(1)=zstore_dummy( ibool(1,NGLLY,1,ispec) ) - xelm(2)=xstore_dummy( ibool(NGLLX,NGLLY,1,ispec) ) - yelm(2)=ystore_dummy( ibool(NGLLX,NGLLY,1,ispec) ) - zelm(2)=zstore_dummy( ibool(NGLLX,NGLLY,1,ispec) ) - xelm(3)=xstore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) - yelm(3)=ystore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) - zelm(3)=zstore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) - xelm(4)=xstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) ) - yelm(4)=ystore_dummy( ibool(1,NGLLY,NGLLZ,ispec) ) - zelm(4)=zstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) ) + xelm(1)=xstore( ibool(1,NGLLY,1,ispec) ) + yelm(1)=ystore( ibool(1,NGLLY,1,ispec) ) + zelm(1)=zstore( ibool(1,NGLLY,1,ispec) ) + xelm(2)=xstore( ibool(NGLLX,NGLLY,1,ispec) ) + yelm(2)=ystore( ibool(NGLLX,NGLLY,1,ispec) ) + zelm(2)=zstore( ibool(NGLLX,NGLLY,1,ispec) ) + xelm(3)=xstore( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) + yelm(3)=ystore( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) + zelm(3)=zstore( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) + xelm(4)=xstore( ibool(1,NGLLY,NGLLZ,ispec) ) + yelm(4)=ystore( ibool(1,NGLLY,NGLLZ,ispec) ) + zelm(4)=zstore( ibool(1,NGLLY,NGLLZ,ispec) ) if (NGNOD2D == 9) then - xelm(5)=xstore_dummy( ibool(MIDX,NGLLY,1,ispec) ) - yelm(5)=ystore_dummy( ibool(MIDX,NGLLY,1,ispec) ) - zelm(5)=zstore_dummy( ibool(MIDX,NGLLY,1,ispec) ) - xelm(6)=xstore_dummy( ibool(NGLLX,NGLLY,MIDZ,ispec) ) - yelm(6)=ystore_dummy( ibool(NGLLX,NGLLY,MIDZ,ispec) ) - zelm(6)=zstore_dummy( ibool(NGLLX,NGLLY,MIDZ,ispec) ) - xelm(7)=xstore_dummy( ibool(MIDX,NGLLY,NGLLZ,ispec) ) - yelm(7)=ystore_dummy( ibool(MIDX,NGLLY,NGLLZ,ispec) ) - zelm(7)=zstore_dummy( ibool(MIDX,NGLLY,NGLLZ,ispec) ) - xelm(8)=xstore_dummy( ibool(1,NGLLY,MIDZ,ispec) ) - yelm(8)=ystore_dummy( ibool(1,NGLLY,MIDZ,ispec) ) - zelm(8)=zstore_dummy( ibool(1,NGLLY,MIDZ,ispec) ) - xelm(9)=xstore_dummy( ibool(MIDX,NGLLY,MIDZ,ispec) ) - yelm(9)=ystore_dummy( ibool(MIDX,NGLLY,MIDZ,ispec) ) - zelm(9)=zstore_dummy( ibool(MIDX,NGLLY,MIDZ,ispec) ) + xelm(5)=xstore( ibool(MIDX,NGLLY,1,ispec) ) + yelm(5)=ystore( ibool(MIDX,NGLLY,1,ispec) ) + zelm(5)=zstore( ibool(MIDX,NGLLY,1,ispec) ) + xelm(6)=xstore( ibool(NGLLX,NGLLY,MIDZ,ispec) ) + yelm(6)=ystore( ibool(NGLLX,NGLLY,MIDZ,ispec) ) + zelm(6)=zstore( ibool(NGLLX,NGLLY,MIDZ,ispec) ) + xelm(7)=xstore( ibool(MIDX,NGLLY,NGLLZ,ispec) ) + yelm(7)=ystore( ibool(MIDX,NGLLY,NGLLZ,ispec) ) + zelm(7)=zstore( ibool(MIDX,NGLLY,NGLLZ,ispec) ) + xelm(8)=xstore( ibool(1,NGLLY,MIDZ,ispec) ) + yelm(8)=ystore( ibool(1,NGLLY,MIDZ,ispec) ) + zelm(8)=zstore( ibool(1,NGLLY,MIDZ,ispec) ) + xelm(9)=xstore( ibool(MIDX,NGLLY,MIDZ,ispec) ) + yelm(9)=ystore( ibool(MIDX,NGLLY,MIDZ,ispec) ) + zelm(9)=zstore( ibool(MIDX,NGLLY,MIDZ,ispec) ) endif call compute_jacobian_2D_face(xelm,yelm,zelm, & @@ -216,74 +216,74 @@ subroutine get_jacobian_boundary_face(nspec, & jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) -! on boundary: bottom + ! on boundary: bottom case (5) - xelm(1)=xstore_dummy( ibool(1,1,1,ispec) ) - yelm(1)=ystore_dummy( ibool(1,1,1,ispec) ) - zelm(1)=zstore_dummy( ibool(1,1,1,ispec) ) - xelm(2)=xstore_dummy( ibool(NGLLX,1,1,ispec) ) - yelm(2)=ystore_dummy( ibool(NGLLX,1,1,ispec) ) - zelm(2)=zstore_dummy( ibool(NGLLX,1,1,ispec) ) - xelm(3)=xstore_dummy( ibool(NGLLX,NGLLY,1,ispec) ) - yelm(3)=ystore_dummy( ibool(NGLLX,NGLLY,1,ispec) ) - zelm(3)=zstore_dummy( ibool(NGLLX,NGLLY,1,ispec) ) - xelm(4)=xstore_dummy( ibool(1,NGLLY,1,ispec) ) - yelm(4)=ystore_dummy( ibool(1,NGLLY,1,ispec) ) - zelm(4)=zstore_dummy( ibool(1,NGLLY,1,ispec) ) + xelm(1)=xstore( ibool(1,1,1,ispec) ) + yelm(1)=ystore( ibool(1,1,1,ispec) ) + zelm(1)=zstore( ibool(1,1,1,ispec) ) + xelm(2)=xstore( ibool(NGLLX,1,1,ispec) ) + yelm(2)=ystore( ibool(NGLLX,1,1,ispec) ) + zelm(2)=zstore( ibool(NGLLX,1,1,ispec) ) + xelm(3)=xstore( ibool(NGLLX,NGLLY,1,ispec) ) + yelm(3)=ystore( ibool(NGLLX,NGLLY,1,ispec) ) + zelm(3)=zstore( ibool(NGLLX,NGLLY,1,ispec) ) + xelm(4)=xstore( ibool(1,NGLLY,1,ispec) ) + yelm(4)=ystore( ibool(1,NGLLY,1,ispec) ) + zelm(4)=zstore( ibool(1,NGLLY,1,ispec) ) if (NGNOD2D == 9) then - xelm(5)=xstore_dummy( ibool(MIDX,1,1,ispec) ) - yelm(5)=ystore_dummy( ibool(MIDX,1,1,ispec) ) - zelm(5)=zstore_dummy( ibool(MIDX,1,1,ispec) ) - xelm(6)=xstore_dummy( ibool(NGLLX,MIDY,1,ispec) ) - yelm(6)=ystore_dummy( ibool(NGLLX,MIDY,1,ispec) ) - zelm(6)=zstore_dummy( ibool(NGLLX,MIDY,1,ispec) ) - xelm(7)=xstore_dummy( ibool(MIDX,NGLLY,1,ispec) ) - yelm(7)=ystore_dummy( ibool(MIDX,NGLLY,1,ispec) ) - zelm(7)=zstore_dummy( ibool(MIDX,NGLLY,1,ispec) ) - xelm(8)=xstore_dummy( ibool(1,MIDY,1,ispec) ) - yelm(8)=ystore_dummy( ibool(1,MIDY,1,ispec) ) - zelm(8)=zstore_dummy( ibool(1,MIDY,1,ispec) ) - xelm(9)=xstore_dummy( ibool(MIDX,MIDY,1,ispec) ) - yelm(9)=ystore_dummy( ibool(MIDX,MIDY,1,ispec) ) - zelm(9)=zstore_dummy( ibool(MIDX,MIDY,1,ispec) ) + xelm(5)=xstore( ibool(MIDX,1,1,ispec) ) + yelm(5)=ystore( ibool(MIDX,1,1,ispec) ) + zelm(5)=zstore( ibool(MIDX,1,1,ispec) ) + xelm(6)=xstore( ibool(NGLLX,MIDY,1,ispec) ) + yelm(6)=ystore( ibool(NGLLX,MIDY,1,ispec) ) + zelm(6)=zstore( ibool(NGLLX,MIDY,1,ispec) ) + xelm(7)=xstore( ibool(MIDX,NGLLY,1,ispec) ) + yelm(7)=ystore( ibool(MIDX,NGLLY,1,ispec) ) + zelm(7)=zstore( ibool(MIDX,NGLLY,1,ispec) ) + xelm(8)=xstore( ibool(1,MIDY,1,ispec) ) + yelm(8)=ystore( ibool(1,MIDY,1,ispec) ) + zelm(8)=zstore( ibool(1,MIDY,1,ispec) ) + xelm(9)=xstore( ibool(MIDX,MIDY,1,ispec) ) + yelm(9)=ystore( ibool(MIDX,MIDY,1,ispec) ) + zelm(9)=zstore( ibool(MIDX,MIDY,1,ispec) ) endif call compute_jacobian_2D_face(xelm,yelm,zelm, & dershape2D_bottom,wgllwgll_xy, & jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) -! on boundary: top + ! on boundary: top case (6) - xelm(1)=xstore_dummy( ibool(1,1,NGLLZ,ispec) ) - yelm(1)=ystore_dummy( ibool(1,1,NGLLZ,ispec) ) - zelm(1)=zstore_dummy( ibool(1,1,NGLLZ,ispec) ) - xelm(2)=xstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) ) - yelm(2)=ystore_dummy( ibool(NGLLX,1,NGLLZ,ispec) ) - zelm(2)=zstore_dummy( ibool(NGLLX,1,NGLLZ,ispec) ) - xelm(3)=xstore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) - yelm(3)=ystore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) - zelm(3)=zstore_dummy( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) - xelm(4)=xstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) ) - yelm(4)=ystore_dummy( ibool(1,NGLLY,NGLLZ,ispec) ) - zelm(4)=zstore_dummy( ibool(1,NGLLY,NGLLZ,ispec) ) + xelm(1)=xstore( ibool(1,1,NGLLZ,ispec) ) + yelm(1)=ystore( ibool(1,1,NGLLZ,ispec) ) + zelm(1)=zstore( ibool(1,1,NGLLZ,ispec) ) + xelm(2)=xstore( ibool(NGLLX,1,NGLLZ,ispec) ) + yelm(2)=ystore( ibool(NGLLX,1,NGLLZ,ispec) ) + zelm(2)=zstore( ibool(NGLLX,1,NGLLZ,ispec) ) + xelm(3)=xstore( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) + yelm(3)=ystore( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) + zelm(3)=zstore( ibool(NGLLX,NGLLY,NGLLZ,ispec) ) + xelm(4)=xstore( ibool(1,NGLLY,NGLLZ,ispec) ) + yelm(4)=ystore( ibool(1,NGLLY,NGLLZ,ispec) ) + zelm(4)=zstore( ibool(1,NGLLY,NGLLZ,ispec) ) if (NGNOD2D == 9) then - xelm(5)=xstore_dummy( ibool(MIDX,1,NGLLZ,ispec) ) - yelm(5)=ystore_dummy( ibool(MIDX,1,NGLLZ,ispec) ) - zelm(5)=zstore_dummy( ibool(MIDX,1,NGLLZ,ispec) ) - xelm(6)=xstore_dummy( ibool(NGLLX,MIDY,NGLLZ,ispec) ) - yelm(6)=ystore_dummy( ibool(NGLLX,MIDY,NGLLZ,ispec) ) - zelm(6)=zstore_dummy( ibool(NGLLX,MIDY,NGLLZ,ispec) ) - xelm(7)=xstore_dummy( ibool(MIDX,NGLLY,NGLLZ,ispec) ) - yelm(7)=ystore_dummy( ibool(MIDX,NGLLY,NGLLZ,ispec) ) - zelm(7)=zstore_dummy( ibool(MIDX,NGLLY,NGLLZ,ispec) ) - xelm(8)=xstore_dummy( ibool(1,MIDY,NGLLZ,ispec) ) - yelm(8)=ystore_dummy( ibool(1,MIDY,NGLLZ,ispec) ) - zelm(8)=zstore_dummy( ibool(1,MIDY,NGLLZ,ispec) ) - xelm(9)=xstore_dummy( ibool(MIDX,MIDY,NGLLZ,ispec) ) - yelm(9)=ystore_dummy( ibool(MIDX,MIDY,NGLLZ,ispec) ) - zelm(9)=zstore_dummy( ibool(MIDX,MIDY,NGLLZ,ispec) ) + xelm(5)=xstore( ibool(MIDX,1,NGLLZ,ispec) ) + yelm(5)=ystore( ibool(MIDX,1,NGLLZ,ispec) ) + zelm(5)=zstore( ibool(MIDX,1,NGLLZ,ispec) ) + xelm(6)=xstore( ibool(NGLLX,MIDY,NGLLZ,ispec) ) + yelm(6)=ystore( ibool(NGLLX,MIDY,NGLLZ,ispec) ) + zelm(6)=zstore( ibool(NGLLX,MIDY,NGLLZ,ispec) ) + xelm(7)=xstore( ibool(MIDX,NGLLY,NGLLZ,ispec) ) + yelm(7)=ystore( ibool(MIDX,NGLLY,NGLLZ,ispec) ) + zelm(7)=zstore( ibool(MIDX,NGLLY,NGLLZ,ispec) ) + xelm(8)=xstore( ibool(1,MIDY,NGLLZ,ispec) ) + yelm(8)=ystore( ibool(1,MIDY,NGLLZ,ispec) ) + zelm(8)=zstore( ibool(1,MIDY,NGLLZ,ispec) ) + xelm(9)=xstore( ibool(MIDX,MIDY,NGLLZ,ispec) ) + yelm(9)=ystore( ibool(MIDX,MIDY,NGLLZ,ispec) ) + zelm(9)=zstore( ibool(MIDX,MIDY,NGLLZ,ispec) ) endif call compute_jacobian_2D_face(xelm,yelm,zelm, & diff --git a/src/shared/write_VTK_data.f90 b/src/shared/write_VTK_data.f90 index 7640837e8..065a17c45 100644 --- a/src/shared/write_VTK_data.f90 +++ b/src/shared/write_VTK_data.f90 @@ -45,7 +45,7 @@ ! maybe we find a better way in future... subroutine write_VTK_data_elem_i(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & elem_flag,prname_file) use constants, only: MAX_STRING_LEN,IOUT_VTK,CUSTOM_REAL,NGLLX,NGLLY,NGLLZ @@ -56,7 +56,7 @@ subroutine write_VTK_data_elem_i(nspec,nglob, & ! global coordinates integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool - real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore ! element flag array integer, dimension(nspec) :: elem_flag @@ -79,7 +79,7 @@ subroutine write_VTK_data_elem_i(nspec,nglob, & write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID' write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float' do i=1,nglob - write(IOUT_VTK,'(3e20.12)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e20.12)') xstore(i),ystore(i),zstore(i) enddo write(IOUT_VTK,*) '' @@ -115,7 +115,7 @@ end subroutine write_VTK_data_elem_i ! routine for saving vtk file holding logical flag on each spectral element subroutine write_VTK_data_elem_l(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & elem_flag,prname_file) use constants, only: MAX_STRING_LEN,IOUT_VTK,CUSTOM_REAL,NGLLX,NGLLY,NGLLZ @@ -126,7 +126,7 @@ subroutine write_VTK_data_elem_l(nspec,nglob, & ! global coordinates integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool - real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore ! element flag array logical, dimension(nspec) :: elem_flag @@ -147,7 +147,7 @@ subroutine write_VTK_data_elem_l(nspec,nglob, & write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID' write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float' do i=1,nglob - write(IOUT_VTK,'(3e20.12)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e20.12)') xstore(i),ystore(i),zstore(i) enddo write(IOUT_VTK,*) '' @@ -188,7 +188,7 @@ end subroutine write_VTK_data_elem_l ! external mesh routine for saving vtk files for CUSTOM_REAL values on all GLL points subroutine write_VTK_data_gll_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & gll_data,prname_file) use constants, only: MAX_STRING_LEN,IOUT_VTK,CUSTOM_REAL,NGLLX,NGLLY,NGLLZ @@ -199,7 +199,7 @@ subroutine write_VTK_data_gll_cr(nspec,nglob, & ! global coordinates integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool - real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore ! GLL data values array real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: gll_data @@ -225,7 +225,7 @@ subroutine write_VTK_data_gll_cr(nspec,nglob, & write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID' write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float' do i=1,nglob - write(IOUT_VTK,'(3e20.12)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e20.12)') xstore(i),ystore(i),zstore(i) enddo write(IOUT_VTK,*) '' @@ -283,7 +283,7 @@ end subroutine write_VTK_data_gll_cr ! external mesh routine for saving vtk files for integer values on all GLL points subroutine write_VTK_data_gll_i(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & gll_data,prname_file) use constants, only: MAX_STRING_LEN,IOUT_VTK,CUSTOM_REAL,NGLLX,NGLLY,NGLLZ @@ -294,7 +294,7 @@ subroutine write_VTK_data_gll_i(nspec,nglob, & ! global coordinates integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool - real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore ! GLL data values array integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: gll_data @@ -320,7 +320,7 @@ subroutine write_VTK_data_gll_i(nspec,nglob, & write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID' write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float' do i=1,nglob - write(IOUT_VTK,'(3e20.12)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e20.12)') xstore(i),ystore(i),zstore(i) enddo write(IOUT_VTK,*) '' @@ -376,7 +376,7 @@ end subroutine write_VTK_data_gll_i ! external mesh routine for saving vtk files for points locations subroutine write_VTK_data_points(nglob, & - xstore_dummy,ystore_dummy,zstore_dummy, & + xstore,ystore,zstore, & points_globalindices,num_points_globalindices, & prname_file) @@ -387,7 +387,7 @@ subroutine write_VTK_data_points(nglob, & integer :: nglob ! global coordinates - real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore ! GLL data values array integer :: num_points_globalindices @@ -418,7 +418,7 @@ subroutine write_VTK_data_points(nglob, & stop 'error vtk points file' endif - write(IOUT_VTK,'(3e20.12)') xstore_dummy(iglob),ystore_dummy(iglob),zstore_dummy(iglob) + write(IOUT_VTK,'(3e20.12)') xstore(iglob),ystore(iglob),zstore(iglob) enddo write(IOUT_VTK,*) '' close(IOUT_VTK) @@ -549,7 +549,7 @@ end subroutine write_VTK_data_points_elem ! external mesh routine for saving vtk files for points locations subroutine write_VTK_data_elem_vectors(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & elem_vector,prname_file) use constants, only: MAX_STRING_LEN,IOUT_VTK,CUSTOM_REAL,NGLLX,NGLLY,NGLLZ @@ -560,7 +560,7 @@ subroutine write_VTK_data_elem_vectors(nspec,nglob, & ! global coordinates integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool - real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore ! element flag array real(kind=CUSTOM_REAL), dimension(3,nspec) :: elem_vector @@ -581,7 +581,7 @@ subroutine write_VTK_data_elem_vectors(nspec,nglob, & write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID' write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float' do i=1,nglob - write(IOUT_VTK,'(3e20.12)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e20.12)') xstore(i),ystore(i),zstore(i) enddo write(IOUT_VTK,*) '' @@ -615,7 +615,7 @@ end subroutine write_VTK_data_elem_vectors ! subroutine write_VTK_data_elem_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & elem_data,prname_file) use constants, only: CUSTOM_REAL,MAX_STRING_LEN,NGLLX,NGLLY,NGLLZ,IOUT_VTK @@ -626,7 +626,7 @@ subroutine write_VTK_data_elem_cr(nspec,nglob, & ! global coordinates integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool - real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore ! element flag array real(kind=CUSTOM_REAL), dimension(nspec) :: elem_data @@ -650,7 +650,7 @@ subroutine write_VTK_data_elem_cr(nspec,nglob, & write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID' write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float' do i = 1,nglob - write(IOUT_VTK,'(3e20.12)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e20.12)') xstore(i),ystore(i),zstore(i) enddo write(IOUT_VTK,*) '' @@ -685,7 +685,7 @@ end subroutine write_VTK_data_elem_cr ! subroutine write_VTU_data_elem_cr_binary(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & elem_data,prname_file) ! saves vtu files in binary format, for CUSTOM_REAL values on all spectral elements @@ -698,7 +698,7 @@ subroutine write_VTU_data_elem_cr_binary(nspec,nglob, & ! global coordinates integer, dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: ibool - real(kind=CUSTOM_REAL), dimension(nglob),intent(in) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL), dimension(nglob),intent(in) :: xstore,ystore,zstore ! element data values array real(kind=CUSTOM_REAL), dimension(nspec),intent(in) :: elem_data @@ -832,7 +832,7 @@ subroutine write_VTU_data_elem_cr_binary(nspec,nglob, & len_bytes = 3 * nglob * SIZE_REAL write(IOUT_VTK) len_bytes do i = 1,nglob - write(IOUT_VTK) real(xstore_dummy(i),kind=4),real(ystore_dummy(i),kind=4),real(zstore_dummy(i),kind=4) + write(IOUT_VTK) real(xstore(i),kind=4),real(ystore(i),kind=4),real(zstore(i),kind=4) enddo ! cells ! connectivity @@ -877,7 +877,7 @@ end subroutine write_VTU_data_elem_cr_binary ! subroutine write_VTK_data_ngnod_elem_i(nspec,nglob,NGNOD, & - xstore_dummy,ystore_dummy,zstore_dummy,elmnts, & + xstore_db,ystore_db,zstore_db,elmnts, & elem_flag,filename) use constants, only: IOUT_VTK,MAX_STRING_LEN @@ -888,7 +888,7 @@ subroutine write_VTK_data_ngnod_elem_i(nspec,nglob,NGNOD, & ! global coordinates integer, dimension(NGNOD,nspec) :: elmnts - double precision, dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + double precision, dimension(nglob) :: xstore_db,ystore_db,zstore_db ! element flag array integer, dimension(nspec) :: elem_flag @@ -910,7 +910,7 @@ subroutine write_VTK_data_ngnod_elem_i(nspec,nglob,NGNOD, & write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID' write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float' do i = 1,nglob - write(IOUT_VTK,'(3e20.12)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e20.12)') xstore_db(i),ystore_db(i),zstore_db(i) enddo write(IOUT_VTK,*) "" @@ -946,7 +946,7 @@ end subroutine write_VTK_data_ngnod_elem_i ! !------------------------------------------------------------------------------------ - subroutine write_VTK_data_elem_cr_meshfem(nspec,nglob,x_dummy,y_dummy,z_dummy,ibool, & + subroutine write_VTK_data_elem_cr_meshfem(nspec,nglob,xstore_db,ystore_db,zstore_db,ibool, & elem_data,filename) ! special routine for meshfem3D with simpler mesh arrays @@ -959,7 +959,7 @@ subroutine write_VTK_data_elem_cr_meshfem(nspec,nglob,x_dummy,y_dummy,z_dummy,ib ! global coordinates integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC) :: ibool - double precision, dimension(nglob) :: x_dummy,y_dummy,z_dummy + double precision, dimension(nglob) :: xstore_db,ystore_db,zstore_db ! element flag array real(kind=CUSTOM_REAL), dimension(nspec) :: elem_data @@ -977,7 +977,7 @@ subroutine write_VTK_data_elem_cr_meshfem(nspec,nglob,x_dummy,y_dummy,z_dummy,ib write(IOUT_VTK,'(a)') 'DATASET UNSTRUCTURED_GRID' write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nglob, ' float' do i = 1,nglob - write(IOUT_VTK,*) sngl(x_dummy(i)),sngl(y_dummy(i)),sngl(z_dummy(i)) + write(IOUT_VTK,*) sngl(xstore_db(i)),sngl(ystore_db(i)),sngl(zstore_db(i)) enddo write(IOUT_VTK,*) '' diff --git a/src/specfem3D/setup_sources_receivers.f90 b/src/specfem3D/setup_sources_receivers.f90 index 4a848875b..c968c2025 100644 --- a/src/specfem3D/setup_sources_receivers.f90 +++ b/src/specfem3D/setup_sources_receivers.f90 @@ -91,6 +91,7 @@ subroutine setup_point_search_arrays() ! prepares midpoints coordinates allocate(xyz_midpoints(NDIM,NSPEC_AB),stat=ier) if (ier /= 0 ) call exit_MPI(myrank,'Error allocating array xyz_midpoints') + xyz_midpoints(:,:) = 0.d0 ! store x/y/z coordinates of center point do ispec = 1,NSPEC_AB @@ -103,6 +104,8 @@ subroutine setup_point_search_arrays() ! define (i,j,k) indices of the control/anchor points allocate(anchor_iax(NGNOD),anchor_iay(NGNOD),anchor_iaz(NGNOD),stat=ier) if (ier /= 0 ) call exit_MPI(myrank,'Error allocating array anchor_i**') + anchor_iax(:) = 0; anchor_iay(:) = 0; anchor_iaz(:) = 0 + call hex_nodes_anchor_ijk_NGLL(NGNOD,anchor_iax,anchor_iay,anchor_iaz,NGLLX,NGLLY,NGLLZ) ! builds search tree @@ -1755,6 +1758,7 @@ subroutine setup_search_kdtree() allocate(kdtree_nodes_location(NDIM,kdtree_num_nodes),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 2096') if (ier /= 0) stop 'Error allocating kdtree_nodes_location arrays' + allocate(kdtree_nodes_index(kdtree_num_nodes),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 2097') if (ier /= 0) stop 'Error allocating kdtree_nodes_index arrays' diff --git a/src/tomography/save_external_bin_m_up.f90 b/src/tomography/save_external_bin_m_up.f90 index 78c41a6f6..e588c4a72 100644 --- a/src/tomography/save_external_bin_m_up.f90 +++ b/src/tomography/save_external_bin_m_up.f90 @@ -37,7 +37,7 @@ subroutine save_external_bin_m_up(nspec,nglob, & rmass,rmass_acoustic,rmass_solid_poroelastic,rmass_fluid_poroelastic, & APPROXIMATE_OCEAN_LOAD,rmass_ocean_load,NGLOB_OCEAN, & ibool, & - xstore_dummy,ystore_dummy,zstore_dummy, & + xstore,ystore,zstore, & abs_boundary_normal,abs_boundary_jacobian2Dw, & abs_boundary_ijk,abs_boundary_ispec, & num_abs_boundary_faces, & @@ -97,7 +97,7 @@ subroutine save_external_bin_m_up(nspec,nglob, & ! mesh coordinates integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool - real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore ! absorbing boundary surface integer :: num_abs_boundary_faces @@ -158,9 +158,9 @@ subroutine save_external_bin_m_up(nspec,nglob, & write(IOUT) ibool - write(IOUT) xstore_dummy - write(IOUT) ystore_dummy - write(IOUT) zstore_dummy + write(IOUT) xstore + write(IOUT) ystore + write(IOUT) zstore write(IOUT) irregular_element_number write(IOUT) xix_regular @@ -333,17 +333,17 @@ subroutine save_external_bin_m_up(nspec,nglob, & ! mesh arrays used for example in combine_vol_data.f90 !--- x coordinate open(unit=IOUT,file=prname(1:len_trim(prname))//'x.bin',status='unknown',form='unformatted') - write(IOUT) xstore_dummy + write(IOUT) xstore close(IOUT) !--- y coordinate open(unit=IOUT,file=prname(1:len_trim(prname))//'y.bin',status='unknown',form='unformatted') - write(IOUT) ystore_dummy + write(IOUT) ystore close(IOUT) !--- z coordinate open(unit=IOUT,file=prname(1:len_trim(prname))//'z.bin',status='unknown',form='unformatted') - write(IOUT) zstore_dummy + write(IOUT) zstore close(IOUT) ! ibool @@ -372,7 +372,7 @@ subroutine save_external_bin_m_up(nspec,nglob, & ! vp values filename = prname(1:len_trim(prname))//'vp_new' call write_VTK_data_gll_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & v_tmp,filename) @@ -393,7 +393,7 @@ subroutine save_external_bin_m_up(nspec,nglob, & ! vs values filename = prname(1:len_trim(prname))//'vs_new' call write_VTK_data_gll_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & v_tmp,filename) ! outputs density model for check @@ -407,7 +407,7 @@ subroutine save_external_bin_m_up(nspec,nglob, & ! density model filename = prname(1:len_trim(prname))//'rho_new' call write_VTK_data_gll_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & v_tmp,filename) @@ -416,7 +416,7 @@ subroutine save_external_bin_m_up(nspec,nglob, & ! saves attenuation flag assigned on each GLL point into a vtk file filename = prname(1:len_trim(prname))//'attenuation' call write_VTK_data_gll_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & qmu_attenuation_store,filename) endif @@ -439,7 +439,7 @@ subroutine save_external_bin_m_up(nspec,nglob, & enddo filename = prname(1:len_trim(prname))//'coupling_acoustic_elastic' call write_VTK_data_points(nglob, & - xstore_dummy,ystore_dummy,zstore_dummy, & + xstore,ystore,zstore, & iglob_tmp,NGLLSQUARE*num_coupling_ac_el_faces, & filename) @@ -457,7 +457,7 @@ subroutine save_external_bin_m_up(nspec,nglob, & enddo filename = prname(1:len_trim(prname))//'acoustic_elastic_flag' call write_VTK_data_elem_i(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & + xstore,ystore,zstore,ibool, & v_tmp_i,filename) endif @@ -465,7 +465,7 @@ subroutine save_external_bin_m_up(nspec,nglob, & ! if (num_interfaces_ext_mesh >= 1) then ! filename = prname(1:len_trim(prname))//'MPI_1_points' ! call write_VTK_data_points(nglob, & - ! xstore_dummy,ystore_dummy,zstore_dummy, & + ! xstore,ystore,zstore, & ! ibool_interfaces_ext_mesh_dummy(1:nibool_interfaces_ext_mesh(1),1), & ! nibool_interfaces_ext_mesh(1), & ! filename) diff --git a/utils/Visualization/Paraview/create_slice_VTK.f90 b/utils/Visualization/Paraview/create_slice_VTK.f90 index 4c58c6ce2..4b73e4041 100644 --- a/utils/Visualization/Paraview/create_slice_VTK.f90 +++ b/utils/Visualization/Paraview/create_slice_VTK.f90 @@ -186,8 +186,8 @@ end program create_slice_VTK ! subroutine write_VTK_data_gll_cr(nspec,nglob, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool, & - gll_data,prname_file) + xstore,ystore,zstore,ibool, & + gll_data,prname_file) ! external mesh routine for saving vtk files for CUSTOM_REAL values on all GLL points @@ -199,7 +199,7 @@ subroutine write_VTK_data_gll_cr(nspec,nglob, & ! global coordinates integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool - real(kind=CUSTOM_REAL), dimension(nglob) :: xstore_dummy,ystore_dummy,zstore_dummy + real(kind=CUSTOM_REAL), dimension(nglob) :: xstore,ystore,zstore ! GLL data values array real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: gll_data @@ -223,28 +223,28 @@ subroutine write_VTK_data_gll_cr(nspec,nglob, & write(IOUT_VTK, '(a,i12,a)') 'POINTS ', nspec*8, ' float' do ispec=1,nspec i = ibool(1,1,1,ispec) - write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e18.6)') xstore(i),ystore(i),zstore(i) i = ibool(NGLLX,1,1,ispec) - write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e18.6)') xstore(i),ystore(i),zstore(i) i = ibool(NGLLX,NGLLY,1,ispec) - write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e18.6)') xstore(i),ystore(i),zstore(i) i = ibool(1,NGLLY,1,ispec) - write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e18.6)') xstore(i),ystore(i),zstore(i) i = ibool(1,1,NGLLZ,ispec) - write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e18.6)') xstore(i),ystore(i),zstore(i) i = ibool(NGLLX,1,NGLLZ,ispec) - write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e18.6)') xstore(i),ystore(i),zstore(i) i = ibool(NGLLX,NGLLY,NGLLZ,ispec) - write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e18.6)') xstore(i),ystore(i),zstore(i) i = ibool(1,NGLLY,NGLLZ,ispec) - write(IOUT_VTK,'(3e18.6)') xstore_dummy(i),ystore_dummy(i),zstore_dummy(i) + write(IOUT_VTK,'(3e18.6)') xstore(i),ystore(i),zstore(i) enddo write(IOUT_VTK,*) ''