From 810b80b80563bab587bf98b928d384c725238a1e Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Mon, 28 Oct 2024 17:17:56 +0100 Subject: [PATCH 1/5] removes unused parameter --- src/generate_databases/setup_mesh_adjacency.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generate_databases/setup_mesh_adjacency.f90 b/src/generate_databases/setup_mesh_adjacency.f90 index 2587cd9d4..932ea4b69 100644 --- a/src/generate_databases/setup_mesh_adjacency.f90 +++ b/src/generate_databases/setup_mesh_adjacency.f90 @@ -31,7 +31,7 @@ subroutine setup_mesh_adjacency() ! setups mesh adjacency array to search element neighbors for point searches use constants, only: myrank, & - NDIM,NGLLX,NGLLY,NGLLZ,MIDX,MIDY,MIDZ,IMAIN,CUSTOM_REAL,MAX_STRING_LEN + NGLLX,NGLLY,NGLLZ,MIDX,MIDY,MIDZ,IMAIN,MAX_STRING_LEN use generate_databases_par, only: NSPEC_AB,NGLOB_AB,ibool,NPROC,prname From d8561f1d91739162a963e76c74e62904c633f623 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Wed, 30 Oct 2024 14:12:12 +0100 Subject: [PATCH 2/5] fixes point search in neighbor elements --- src/shared/recompute_jacobian.f90 | 63 ++++++++------ src/specfem3D/locate_point.f90 | 132 ++++++++++++++++++++--------- src/specfem3D/locate_receivers.f90 | 36 +++++++- 3 files changed, 166 insertions(+), 65 deletions(-) diff --git a/src/shared/recompute_jacobian.f90 b/src/shared/recompute_jacobian.f90 index 7ea10e8ba..1aa9c26e1 100644 --- a/src/shared/recompute_jacobian.f90 +++ b/src/shared/recompute_jacobian.f90 @@ -30,7 +30,7 @@ subroutine recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, & xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,NGNOD) - use constants, only: NDIM,ZERO + use constants, only: NDIM,ZERO,myrank implicit none @@ -56,7 +56,7 @@ subroutine recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, & ! recompute jacobian for any (xi,eta,gamma) point, not necessarily a GLL point -! check that the parameter file is correct + ! check that the parameter file is correct if (NGNOD /= 8 .and. NGNOD /= 27) stop 'elements should have 8 or 27 control nodes' if (NGNOD == 8) then @@ -71,51 +71,66 @@ subroutine recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, & x = ZERO y = ZERO z = ZERO + xxi = ZERO xeta = ZERO xgamma = ZERO + yxi = ZERO yeta = ZERO ygamma = ZERO + zxi = ZERO zeta = ZERO zgamma = ZERO do ia = 1,NGNOD - x = x+shape3D(ia)*xelm(ia) - y = y+shape3D(ia)*yelm(ia) - z = z+shape3D(ia)*zelm(ia) + x = x + shape3D(ia) * xelm(ia) + y = y + shape3D(ia) * yelm(ia) + z = z + shape3D(ia) * zelm(ia) ! Jacobian matrix elements ! dx/dxi, dx/deta, etc. - xxi = xxi+dershape3D(1,ia)*xelm(ia) - xeta = xeta+dershape3D(2,ia)*xelm(ia) - xgamma = xgamma+dershape3D(3,ia)*xelm(ia) - yxi = yxi+dershape3D(1,ia)*yelm(ia) - yeta = yeta+dershape3D(2,ia)*yelm(ia) - ygamma = ygamma+dershape3D(3,ia)*yelm(ia) - zxi = zxi+dershape3D(1,ia)*zelm(ia) - zeta = zeta+dershape3D(2,ia)*zelm(ia) - zgamma = zgamma+dershape3D(3,ia)*zelm(ia) + xxi = xxi + dershape3D(1,ia)*xelm(ia) + xeta = xeta + dershape3D(2,ia)*xelm(ia) + xgamma = xgamma + dershape3D(3,ia)*xelm(ia) + + yxi = yxi + dershape3D(1,ia)*yelm(ia) + yeta = yeta + dershape3D(2,ia)*yelm(ia) + ygamma = ygamma + dershape3D(3,ia)*yelm(ia) + + zxi = zxi + dershape3D(1,ia)*zelm(ia) + zeta = zeta + dershape3D(2,ia)*zelm(ia) + zgamma = zgamma + dershape3D(3,ia)*zelm(ia) enddo ! Jacobian determinant jacobian = xxi*(yeta*zgamma-ygamma*zeta) - xeta*(yxi*zgamma-ygamma*zxi) + xgamma*(yxi*zeta-yeta*zxi) - if (jacobian <= ZERO) stop '3D Jacobian undefined' + if (jacobian <= ZERO) then + ! info + print *,'Error Jacobian rank:',myrank,'has invalid Jacobian ',jacobian + print *,' input xi/eta/gamma: ',xi,eta,gamma + print *,' Jacobian: ',jacobian,'xxi',xxi,xeta,xgamma,'yxi',yxi,yeta,ygamma,'zxi',zxi,zeta,zgamma + print *,' location x/y/z: ',x,y,z + + stop '3D Jacobian undefined' + endif ! inverse Jacobian matrix ! ! invert the relation (Fletcher p. 50 vol. 2) - xix = (yeta*zgamma-ygamma*zeta)/jacobian - xiy = (xgamma*zeta-xeta*zgamma)/jacobian - xiz = (xeta*ygamma-xgamma*yeta)/jacobian - etax = (ygamma*zxi-yxi*zgamma)/jacobian - etay = (xxi*zgamma-xgamma*zxi)/jacobian - etaz = (xgamma*yxi-xxi*ygamma)/jacobian - gammax = (yxi*zeta-yeta*zxi)/jacobian - gammay = (xeta*zxi-xxi*zeta)/jacobian - gammaz = (xxi*yeta-xeta*yxi)/jacobian + xix = (yeta*zgamma - ygamma*zeta) / jacobian + xiy = (xgamma*zeta - xeta*zgamma) / jacobian + xiz = (xeta*ygamma - xgamma*yeta) / jacobian + + etax = (ygamma*zxi - yxi*zgamma) / jacobian + etay = (xxi*zgamma - xgamma*zxi) / jacobian + etaz = (xgamma*yxi - xxi*ygamma) / jacobian + + gammax = (yxi*zeta - yeta*zxi) / jacobian + gammay = (xeta*zxi - xxi*zeta) / jacobian + gammaz = (xxi*yeta - xeta*yxi) / jacobian end subroutine recompute_jacobian diff --git a/src/specfem3D/locate_point.f90 b/src/specfem3D/locate_point.f90 index 95156aca7..9beeae218 100644 --- a/src/specfem3D/locate_point.f90 +++ b/src/specfem3D/locate_point.f90 @@ -32,7 +32,8 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & POINT_CAN_BE_BURIED, elemsize_max_glob, & ispec_selected, xi_found, eta_found, gamma_found, & - x_found, y_found, z_found, domain, nu_point, final_distance_squared) + x_found, y_found, z_found, & + domain, nu_point, final_distance_squared) use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM, & IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,IDOMAIN_POROELASTIC, & @@ -69,6 +70,7 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & double precision :: distmin_squared,dist_squared double precision :: x,y,z double precision :: xi,eta,gamma + double precision :: dx,dy,dz integer :: ix_initial_guess, iy_initial_guess, iz_initial_guess @@ -152,6 +154,8 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & y_found = -HUGEVAL z_found = -HUGEVAL + nu_point(:,:) = 0.d0 + ! store final distance squared between asked and found final_distance_squared = HUGEVAL @@ -215,9 +219,10 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & endif ! distance to element midpoint - dist_squared = (x_target - xyz_midpoints(1,ispec))*(x_target - xyz_midpoints(1,ispec)) & - + (y_target - xyz_midpoints(2,ispec))*(y_target - xyz_midpoints(2,ispec)) & - + (z_target - xyz_midpoints(3,ispec))*(z_target - xyz_midpoints(3,ispec)) + dx = x_target - xyz_midpoints(1,ispec) + dy = y_target - xyz_midpoints(2,ispec) + dz = z_target - xyz_midpoints(3,ispec) + dist_squared = dx*dx + dy*dy + dz*dz ! we compare squared distances instead of distances to speed up the comparison by avoiding computing a square root if (dist_squared > maximal_elem_size_squared) cycle ! exclude elements that are too far from target @@ -239,9 +244,10 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & z = dble(zstore(iglob)) ! distance to GLL point - dist_squared = (x_target - x)*(x_target - x) & - + (y_target - y)*(y_target - y) & - + (z_target - z)*(z_target - z) + dx = x_target - x + dy = y_target - y + dz = z_target - z + dist_squared = dx*dx + dy*dy + dz*dz if (dist_squared < distmin_squared) then distmin_squared = dist_squared @@ -317,7 +323,10 @@ subroutine locate_point_in_mesh(x_target, y_target, z_target, & ispec_selected,ix_initial_guess,iy_initial_guess,iz_initial_guess) ! compute distance squared between asked and found - dist_squared = (x_target-x)*(x_target-x) + (y_target-y)*(y_target-y) + (z_target-z)*(z_target-z) + dx = x_target - x + dy = y_target - y + dz = z_target - z + dist_squared = dx*dx + dy*dy + dz*dz ! initial point estimate xi_found = xi @@ -426,8 +435,6 @@ subroutine find_local_coordinates(x_target,y_target,z_target,xi,eta,gamma,x,y,z, use specfem_par, only: xstore,ystore,zstore,ibool,NGNOD,anchor_iax,anchor_iay,anchor_iaz, & xigll,yigll,zigll - !use constants, only: myrank - implicit none double precision,intent(in) :: x_target,y_target,z_target @@ -531,19 +538,22 @@ subroutine find_local_coordinates(x_target,y_target,z_target,xi,eta,gamma,x,y,z, ! ! note: we are more conservative here than in the 3D_GLOBE version, since the mesh can have voids or complex shapes, ! and we want the point location not to be too far off the mesh volume - if (xi > 1.01d0) xi = 1.01d0 - if (xi < -1.01d0) xi = -1.01d0 - if (eta > 1.01d0) eta = 1.01d0 - if (eta < -1.01d0) eta = -1.01d0 - if (gamma > 1.01d0) gamma = 1.01d0 - if (gamma < -1.01d0) gamma = -1.01d0 + !if (xi > 1.01d0) xi = 1.01d0 + !if (xi <-1.01d0) xi = -1.01d0 + if (abs(xi) > 1.01d0) xi = sign(1.01d0,xi) + !if (eta > 1.01d0) eta = 1.01d0 + !if (eta <-1.01d0) eta = -1.01d0 + if (abs(eta) > 1.01d0) eta = sign(1.01d0,eta) + !if (gamma > 1.01d0) gamma = 1.01d0 + !if (gamma <-1.01d0) gamma = -1.01d0 + if (abs(gamma) > 1.01d0) gamma = sign(1.01d0,gamma) ! end of non linear iterations enddo ! position refinements ! check distance and add more refinement if too far off - if (d_min_sq > 1.0) then + if (d_min_sq > 1.d0) then ! only refine if still within element if (abs(xi) <= 1.d0 .and. abs(eta) <= 1.d0 .and. abs(gamma) <= 1.d0) then ! moves point slightly off current position @@ -596,20 +606,26 @@ subroutine find_local_coordinates(x_target,y_target,z_target,xi,eta,gamma,x,y,z, eta = eta + deta gamma = gamma + dgamma - if (xi > 1.01d0) xi = 1.01d0 - if (xi < -1.01d0) xi = -1.01d0 - if (eta > 1.01d0) eta = 1.01d0 - if (eta < -1.01d0) eta = -1.01d0 - if (gamma > 1.01d0) gamma = 1.01d0 - if (gamma < -1.01d0) gamma = -1.01d0 + !if (xi > 1.01d0) xi = 1.01d0 + !if (xi <-1.01d0) xi = -1.01d0 + if (abs(xi) > 1.01d0) xi = sign(1.01d0,xi) + !if (eta > 1.01d0) eta = 1.01d0 + !if (eta <-1.01d0) eta = -1.01d0 + if (abs(eta) > 1.01d0) eta = sign(1.01d0,eta) + !if (gamma > 1.01d0) gamma = 1.01d0 + !if (gamma <-1.01d0) gamma = -1.01d0 + if (abs(gamma) > 1.01d0) gamma = sign(1.01d0,gamma) ! stop if point moves outside of that element - if (xi >= 1.01d0) exit - if (xi <= -1.01d0) exit - if (eta >= 1.01d0) exit - if (eta <= -1.01d0) exit - if (gamma >= 1.01d0) exit - if (gamma <= -1.01d0) exit + !if (xi >= 1.01d0) exit + !if (xi <= -1.01d0) exit + !if (eta >= 1.01d0) exit + !if (eta <= -1.01d0) exit + !if (gamma >= 1.01d0) exit + !if (gamma <= -1.01d0) exit + if (abs(xi) >= 1.01d0) exit + if (abs(eta) >= 1.01d0) exit + if (abs(gamma) >= 1.01d0) exit ! stop criteria with accuracy below 1.d-10 if (d_min_sq < 1.d-10) exit @@ -644,6 +660,10 @@ subroutine find_best_neighbor(x_target,y_target,z_target, & use specfem_par_acoustic, only: ispec_is_acoustic use specfem_par_elastic, only: ispec_is_elastic + !debug + use constants, only: myrank + use specfem_par, only: NSPEC_AB + implicit none double precision, intent(in) :: x_target, y_target, z_target @@ -657,11 +677,12 @@ subroutine find_best_neighbor(x_target,y_target,z_target, & logical, intent(in) :: POINT_CAN_BE_BURIED ! locals - integer :: ispec,i,j,k,iglob + integer :: ispec_ref,ispec,i,j,k,iglob ! location search double precision :: final_distance_squared_this_element double precision :: x,y,z double precision :: xi,eta,gamma + double precision :: dx,dy,dz ! neighbor elements integer :: num_neighbors @@ -670,17 +691,46 @@ subroutine find_best_neighbor(x_target,y_target,z_target, & double precision :: distmin_squared_guess,dist_squared logical :: is_better_location + ! for close points + double precision, parameter :: TOL_POINT_DISTANCE = 1.d-10 + ! for close points relative to each other + double precision, parameter :: TOL_RELATIVE_POINT_DISTANCE = 1.d-15 + ! note: find_local_coordinates(..) gets first called in any case in the calling routine locate point for an initial guess, ! before running this neighbor point search. ! we can thus skip the reference element (ispec_selected) and only loop over neighbors. + ! reference element is the initial selection + ! we will try to find a better position in these neighbor elements and update the `ispec_selected` selection accordingly + ispec_ref = ispec_selected + ! loops over neighboring elements - num_neighbors = neighbors_xadj(ispec_selected+1) - neighbors_xadj(ispec_selected) + num_neighbors = neighbors_xadj(ispec_ref+1) - neighbors_xadj(ispec_ref) do ii = 1,num_neighbors - ! get neighbor - ientry = neighbors_xadj(ispec_selected) + ii + ! gets neighbor index + ientry = neighbors_xadj(ispec_ref) + ii + + ! checks entry + if (ientry < 1 .or. ientry > size(neighbors_adjncy)) then + print *,'Error: neighbor entry not valid for referenc element ',ispec_ref + print *,' num_neighbors : ',num_neighbors + print *,' neighbor entry : ',ientry,'out of',size(neighbors_adjncy) + print *,' neighbor_xadj : ',neighbors_xadj(ispec_ref),neighbors_xadj(ispec_ref+1) + print *,' neighbor_adjncy: ',neighbors_adjncy(neighbors_xadj(ispec_ref):neighbors_xadj(ispec_ref)+num_neighbors) + call exit_MPI(myrank,'Invalid neighbors_xadj array') + endif + + ! gets element index ispec = neighbors_adjncy(ientry) + ! checks element index + if (ispec < 1 .or. ispec > NSPEC_AB) then + print *,'Error: neighbor not valid for reference element ',ispec_ref + print *,' total neighbors: ',num_neighbors + print *,' neighbor entry: ',ientry,'out of',size(neighbors_adjncy),' has element ispec = ',ispec + call exit_MPI(myrank,'Invalid neighbors_adjncy array') + endif + ! only loop through surface elements if point location must stick to surface if (.not. POINT_CAN_BE_BURIED) then if (.not. ispec_is_surface_external_mesh(ispec)) cycle @@ -713,9 +763,10 @@ subroutine find_best_neighbor(x_target,y_target,z_target, & z = dble(zstore(iglob)) ! distance to GLL point - dist_squared = (x_target - x)*(x_target - x) & - + (y_target - y)*(y_target - y) & - + (z_target - z)*(z_target - z) + dx = x_target - x + dy = y_target - y + dz = z_target - z + dist_squared = dx*dx + dy*dy + dz*dz ! keep this point if it is closer to the receiver ! we compare squared distances instead of distances themselves to significantly speed up calculations @@ -734,13 +785,16 @@ subroutine find_best_neighbor(x_target,y_target,z_target, & ispec,ix_initial_guess,iy_initial_guess,iz_initial_guess) ! compute final distance squared between asked and found - final_distance_squared_this_element = (x_target-x)*(x_target-x) + (y_target-y)*(y_target-y) + (z_target-z)*(z_target-z) + dx = x_target - x + dy = y_target - y + dz = z_target - z + final_distance_squared_this_element = dx*dx + dy*dy + dz*dz ! flag to determine if taking new position is_better_location = .false. ! if we found a close point - if (final_distance_squared_this_element < 1.d-10) then + if (final_distance_squared_this_element < TOL_POINT_DISTANCE) then ! checks if position within element if (abs(xi) <= 1.d0 .and. abs(eta) <= 1.d0 .and. abs(gamma) <= 1.d0) then ! takes position if old position has xi/eta/gamma outside of element @@ -769,7 +823,7 @@ subroutine find_best_neighbor(x_target,y_target,z_target, & ! has been located slightly away from the element. thus, we try to check if the previous location ! was accurate enough with coordinates still within an element. is_better_location = .true. - if (abs(final_distance_squared_this_element - final_distance_squared) < 1.d-15) then + if (abs(final_distance_squared_this_element - final_distance_squared) < TOL_RELATIVE_POINT_DISTANCE) then if (abs(xi) > 1.d0 .or. abs(eta) > 1.d0 .or. abs(gamma) > 1.d0) then if (abs(xi_found) > 1.d0 .or. abs(eta_found) > 1.d0 .or. abs(gamma_found) > 1.d0) then ! old position is also out of element, let's take the new one diff --git a/src/specfem3D/locate_receivers.f90 b/src/specfem3D/locate_receivers.f90 index d8b4cb7b1..b3301cbdf 100644 --- a/src/specfem3D/locate_receivers.f90 +++ b/src/specfem3D/locate_receivers.f90 @@ -283,6 +283,37 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) call bcast_all_dp(nu_rec,NDIM*NDIM*nrec) call bcast_all_dp(final_distance,nrec) + ! checks if valid receiver element + if (myrank == 0) then + ! locate point might return a zero ispec if point outside/above mesh + do irec = 1,nrec + ispec = ispec_selected_rec(irec) + ! checks if valid + if (ispec < 1 .or. ispec > NSPEC_AB) then + ! invalid element + print *,'Error locating station # ',irec,' ',trim(network_name(irec)),' ',trim(station_name(irec)) + if (SUPPRESS_UTM_PROJECTION) then + print *,' original x: ',sngl(stutm_x(irec)) + print *,' original y: ',sngl(stutm_y(irec)) + else + print *,' original UTM x: ',sngl(stutm_x(irec)) + print *,' original UTM y: ',sngl(stutm_y(irec)) + endif + if (USE_SOURCES_RECEIVERS_Z) then + print *,' original z: ',sngl(stbur(irec)) + else + print *,' original depth: ',sngl(stbur(irec)),' m' + endif + print *,' only found invalid element: slice ',islice_selected_rec(irec) + print *,' ispec ',ispec_selected_rec(irec) + print *,' domain ',idomain(irec) + print *,' final_distance: ',final_distance(irec) + call exit_MPI(myrank,'Error locating receiver') + endif + enddo + endif + call synchronize_all() + ! warning if receiver in C-PML region allocate(is_CPML_rec(nrec),stat=ier) if (ier /= 0) call exit_MPI(myrank,'Error allocating is_CPML_rec array') @@ -298,6 +329,7 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) do irec = 1,nrec if (islice_selected_rec(irec) == myrank) then ispec = ispec_selected_rec(irec) + ! check if PML element if (is_CPML(ispec)) then is_CPML_rec(irec) = .true. ! debug @@ -315,8 +347,8 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) ! checks stations location if (final_distance(irec) == HUGEVAL) then - write(IMAIN,*) 'error locating station # ',irec,' ',trim(network_name(irec)),' ',trim(station_name(irec)) - call exit_MPI(myrank,'error locating receiver') + write(IMAIN,*) 'Error locating station # ',irec,' ',trim(network_name(irec)),' ',trim(station_name(irec)) + call exit_MPI(myrank,'Error locating receiver') endif ! limits user output if too many receivers From 49823106c93a45f0a77293cf0c2c508413046896 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Wed, 30 Oct 2024 14:13:26 +0100 Subject: [PATCH 3/5] prefers slices with point locations inside elements --- src/specfem3D/locate_MPI_slice.f90 | 137 ++++++++++++++++++++++++----- 1 file changed, 116 insertions(+), 21 deletions(-) diff --git a/src/specfem3D/locate_MPI_slice.f90 b/src/specfem3D/locate_MPI_slice.f90 index 5305e7f49..92d708a3f 100644 --- a/src/specfem3D/locate_MPI_slice.f90 +++ b/src/specfem3D/locate_MPI_slice.f90 @@ -69,6 +69,17 @@ subroutine locate_MPI_slice(npoints_subset,ipoin_already_done, & double precision, dimension(npoints_subset,0:NPROC-1) :: final_distance_all double precision, dimension(NDIM,NDIM,npoints_subset,0:NPROC-1) :: nu_all + double precision :: xi,eta,gamma + + ! point inside element + logical :: found_point_inside + + ! prefer slice with point location inside element rather than just selecting slice with closest location + logical, parameter :: DO_PREFER_INSIDE_ELEMENT = .true. + + !debug + !integer,dimension(1) :: iproc_min + ! initializes with dummy values ispec_selected_all(:,:) = -1 idomain_all(:,:) = -1000 @@ -98,33 +109,117 @@ subroutine locate_MPI_slice(npoints_subset,ipoin_already_done, & ! find the slice and element to put the source if (myrank == 0) then - ! loops over subset - do ipoin_in_this_subset = 1,npoints_subset + ! we prefer slices with point locations inside an element over slices with a closer point distance + ! but a point location outside the (best) element + if (DO_PREFER_INSIDE_ELEMENT) then + + ! loops over subset + do ipoin_in_this_subset = 1,npoints_subset + ! mapping from station/source number in current subset to real station/source number in all the subsets + ipoin = ipoin_in_this_subset + ipoin_already_done + + ! checks first if we have a close point inside an element + found_point_inside = .false. + distmin = HUGEVAL + do iproc = 0,NPROC-1 + ! point position + xi = xi_all(ipoin_in_this_subset,iproc) + eta = eta_all(ipoin_in_this_subset,iproc) + gamma = gamma_all(ipoin_in_this_subset,iproc) + + ! points inside an element have xi/eta/gamma <= 1.0 + if (abs(xi) <= 1.d0 .and. abs(eta) <= 1.d0 .and. abs(gamma) <= 1.d0) then + ! takes point if closer + if (final_distance_all(ipoin_in_this_subset,iproc) < distmin) then + found_point_inside = .true. + + xi_point(ipoin) = xi + eta_point(ipoin) = eta + gamma_point(ipoin) = gamma + + distmin = final_distance_all(ipoin_in_this_subset,iproc) + + islice_selected(ipoin) = iproc + ispec_selected(ipoin) = ispec_selected_all(ipoin_in_this_subset,iproc) + idomain(ipoin) = idomain_all(ipoin_in_this_subset,iproc) + + nu_point(:,:,ipoin) = nu_all(:,:,ipoin_in_this_subset,iproc) + + x_found(ipoin) = x_found_all(ipoin_in_this_subset,iproc) + y_found(ipoin) = y_found_all(ipoin_in_this_subset,iproc) + z_found(ipoin) = z_found_all(ipoin_in_this_subset,iproc) + endif + endif + enddo + + ! if we haven't found a close point inside an element, then look for just the closest possible + if (.not. found_point_inside) then + do iproc = 0,NPROC-1 + if (final_distance_all(ipoin_in_this_subset,iproc) < distmin) then + distmin = final_distance_all(ipoin_in_this_subset,iproc) + + islice_selected(ipoin) = iproc + ispec_selected(ipoin) = ispec_selected_all(ipoin_in_this_subset,iproc) + idomain(ipoin) = idomain_all(ipoin_in_this_subset,iproc) + + xi_point(ipoin) = xi_all(ipoin_in_this_subset,iproc) + eta_point(ipoin) = eta_all(ipoin_in_this_subset,iproc) + gamma_point(ipoin) = gamma_all(ipoin_in_this_subset,iproc) + + nu_point(:,:,ipoin) = nu_all(:,:,ipoin_in_this_subset,iproc) + + x_found(ipoin) = x_found_all(ipoin_in_this_subset,iproc) + y_found(ipoin) = y_found_all(ipoin_in_this_subset,iproc) + z_found(ipoin) = z_found_all(ipoin_in_this_subset,iproc) + endif + enddo + endif - ! mapping from station/source number in current subset to real station/source number in all the subsets - ipoin = ipoin_in_this_subset + ipoin_already_done + final_distance(ipoin) = distmin + + !debug + !iproc_min = minloc(final_distance_all(ipoin_in_this_subset,:)) - 1 ! -1 to start procs at 0 + !print *,'debug: locate MPI point',ipoin,ipoin_in_this_subset,'dist',final_distance_all(ipoin_in_this_subset,:), & + ! 'iproc min',iproc_min(1), & + ! 'xi min',xi_all(ipoin_in_this_subset,iproc_min(1)), & + ! eta_all(ipoin_in_this_subset,iproc_min(1)), & + ! gamma_all(ipoin_in_this_subset,iproc_min(1)), & + ! 'iproc found',islice_selected(ipoin),'dist found',final_distance(ipoin), & + ! 'xi found',xi_point(ipoin),eta_point(ipoin),gamma_point(ipoin) + enddo - distmin = HUGEVAL - do iproc = 0,NPROC-1 - if (final_distance_all(ipoin_in_this_subset,iproc) < distmin) then - distmin = final_distance_all(ipoin_in_this_subset,iproc) + else + ! old version takes closest point - islice_selected(ipoin) = iproc - ispec_selected(ipoin) = ispec_selected_all(ipoin_in_this_subset,iproc) - idomain(ipoin) = idomain_all(ipoin_in_this_subset,iproc) + ! loops over subset + do ipoin_in_this_subset = 1,npoints_subset - xi_point(ipoin) = xi_all(ipoin_in_this_subset,iproc) - eta_point(ipoin) = eta_all(ipoin_in_this_subset,iproc) - gamma_point(ipoin) = gamma_all(ipoin_in_this_subset,iproc) - nu_point(:,:,ipoin) = nu_all(:,:,ipoin_in_this_subset,iproc) + ! mapping from station/source number in current subset to real station/source number in all the subsets + ipoin = ipoin_in_this_subset + ipoin_already_done - x_found(ipoin) = x_found_all(ipoin_in_this_subset,iproc) - y_found(ipoin) = y_found_all(ipoin_in_this_subset,iproc) - z_found(ipoin) = z_found_all(ipoin_in_this_subset,iproc) - endif + distmin = HUGEVAL + do iproc = 0,NPROC-1 + if (final_distance_all(ipoin_in_this_subset,iproc) < distmin) then + distmin = final_distance_all(ipoin_in_this_subset,iproc) + + islice_selected(ipoin) = iproc + ispec_selected(ipoin) = ispec_selected_all(ipoin_in_this_subset,iproc) + idomain(ipoin) = idomain_all(ipoin_in_this_subset,iproc) + + xi_point(ipoin) = xi_all(ipoin_in_this_subset,iproc) + eta_point(ipoin) = eta_all(ipoin_in_this_subset,iproc) + gamma_point(ipoin) = gamma_all(ipoin_in_this_subset,iproc) + nu_point(:,:,ipoin) = nu_all(:,:,ipoin_in_this_subset,iproc) + + x_found(ipoin) = x_found_all(ipoin_in_this_subset,iproc) + y_found(ipoin) = y_found_all(ipoin_in_this_subset,iproc) + z_found(ipoin) = z_found_all(ipoin_in_this_subset,iproc) + endif + enddo + final_distance(ipoin) = distmin enddo - final_distance(ipoin) = distmin - enddo + + endif endif ! end of section executed by main process only From 394063d8c6a6e88107a07692cc8ec17d1700fe43 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Wed, 30 Oct 2024 14:16:15 +0100 Subject: [PATCH 4/5] explicitly specifies size(**,kind=4) in function arguments (to avoid compiler problems where size() returns either integer(kind=4) or integer(kind=8) depending on flags like -mcmedium) --- .../model_external_values.F90 | 2 +- .../model_salton_trough.f90 | 2 +- src/generate_databases/model_tomography.f90 | 2 +- .../setup_mesh_adjacency.f90 | 2 +- src/meshfem3D/determine_cavity.f90 | 4 +- src/specfem3D/iterate_time.F90 | 4 +- src/specfem3D/prepare_attenuation.f90 | 12 +- src/specfem3D/prepare_gpu.f90 | 4 +- src/specfem3D/read_forward_arrays.f90 | 6 +- src/specfem3D/read_mesh_databases.F90 | 253 ++++++++--------- src/specfem3D/read_mesh_databases_hdf5.F90 | 255 +++++++++--------- src/specfem3D/save_forward_arrays.f90 | 2 +- src/specfem3D/vtk_window.F90 | 4 +- src/specfem3D/write_movie_output_HDF5.F90 | 4 +- 14 files changed, 281 insertions(+), 275 deletions(-) diff --git a/src/generate_databases/model_external_values.F90 b/src/generate_databases/model_external_values.F90 index e50cb312d..9930f0cda 100644 --- a/src/generate_databases/model_external_values.F90 +++ b/src/generate_databases/model_external_values.F90 @@ -94,7 +94,7 @@ subroutine model_external_broadcast() if (myrank == 0) call read_external_model() ! broadcast the information read on the main to the nodes - call bcast_all_dp(MEXT_V%dvs, size(MEXT_V%dvs)) + call bcast_all_dp(MEXT_V%dvs, size(MEXT_V%dvs,kind=4)) end subroutine model_external_broadcast diff --git a/src/generate_databases/model_salton_trough.f90 b/src/generate_databases/model_salton_trough.f90 index 1080691aa..78f294656 100644 --- a/src/generate_databases/model_salton_trough.f90 +++ b/src/generate_databases/model_salton_trough.f90 @@ -71,7 +71,7 @@ subroutine model_salton_trough_broadcast() if (myrank == 0) call read_salton_sea_model() ! broadcast the information read on the main to the nodes - call bcast_all_r(vp_array, size(vp_array)) + call bcast_all_r(vp_array, size(vp_array,kind=4)) end subroutine model_salton_trough_broadcast diff --git a/src/generate_databases/model_tomography.f90 b/src/generate_databases/model_tomography.f90 index 9b0c2059f..f5d4a0cdf 100644 --- a/src/generate_databases/model_tomography.f90 +++ b/src/generate_databases/model_tomography.f90 @@ -121,7 +121,7 @@ subroutine model_tomography_broadcast(myrank) ! allocate( vp_tomography(1:nrecord) ,stat=ier) ! if (ier /= 0) stop 'error allocating array vp_tomography' !endif - !call bcast_all_cr(vp_tomography,size(vp_tomography)) + !call bcast_all_cr(vp_tomography,size(vp_tomography,kind=4)) ! synchronizes processes call synchronize_all() diff --git a/src/generate_databases/setup_mesh_adjacency.f90 b/src/generate_databases/setup_mesh_adjacency.f90 index 932ea4b69..23d2a1ecd 100644 --- a/src/generate_databases/setup_mesh_adjacency.f90 +++ b/src/generate_databases/setup_mesh_adjacency.f90 @@ -377,7 +377,7 @@ subroutine setup_mesh_adjacency() endif ! check if element has neighbors - ! note: in case of a fault in this slice (splitting nodes) and/or scotch paritioning + ! note: in case of a fault in this slice (splitting nodes) and/or scotch partitioning ! it can happen that an element has no neighbors if (NPROC == 1 .and. (.not. ANY_FAULT_IN_THIS_PROC)) then ! checks if neighbors were found diff --git a/src/meshfem3D/determine_cavity.f90 b/src/meshfem3D/determine_cavity.f90 index 4127751d1..4871109ba 100644 --- a/src/meshfem3D/determine_cavity.f90 +++ b/src/meshfem3D/determine_cavity.f90 @@ -366,13 +366,13 @@ subroutine cmm_determine_cavity(nglob) allocate(tmp_all(1,1),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 1330') endif - call sum_all_1Darray_dp(cavity_boundary,tmp_all,size(cavity_boundary)) + call sum_all_1Darray_dp(cavity_boundary,tmp_all,size(cavity_boundary,kind=4)) if (myrank == 0) then cavity_boundary(:,:) = tmp_all(:,:) endif deallocate(tmp_all) ! broadcasts to all others - call bcast_all_dp(cavity_boundary,size(cavity_boundary)) + call bcast_all_dp(cavity_boundary,size(cavity_boundary,kind=4)) !print *,'cavity boundary after:',myrank,'array:',cavity_boundary(:,:) diff --git a/src/specfem3D/iterate_time.F90 b/src/specfem3D/iterate_time.F90 index 5c413af04..05943a71c 100644 --- a/src/specfem3D/iterate_time.F90 +++ b/src/specfem3D/iterate_time.F90 @@ -418,9 +418,9 @@ subroutine it_transfer_from_GPU() if (ATTENUATION) then call transfer_rmemory_from_device(Mesh_pointer,R_xx,R_yy,R_xy,R_xz,R_yz, & - R_trace,size(R_xx)) + R_trace,size(R_xx,kind=4)) call transfer_strain_from_device(Mesh_pointer,epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz, & - epsilondev_trace,size(epsilondev_xx)) + epsilondev_trace,size(epsilondev_xx,kind=4)) endif endif diff --git a/src/specfem3D/prepare_attenuation.f90 b/src/specfem3D/prepare_attenuation.f90 index 0f2056afd..afdc85c1e 100644 --- a/src/specfem3D/prepare_attenuation.f90 +++ b/src/specfem3D/prepare_attenuation.f90 @@ -114,12 +114,12 @@ subroutine prepare_attenuation() ! broadcasts call bcast_all_i_for_database(ispec, 1) - if (size(factor_common) > 0) & - call bcast_all_cr_for_database(factor_common(1,1,1,1,1), size(factor_common)) - if (size(scale_factor) > 0) & - call bcast_all_cr_for_database(scale_factor(1,1,1,1), size(scale_factor)) - call bcast_all_cr_for_database(factor_common_kappa(1,1,1,1,1), size(factor_common_kappa)) - call bcast_all_cr_for_database(scale_factor_kappa(1,1,1,1), size(scale_factor_kappa)) + if (size(factor_common,kind=4) > 0) & + call bcast_all_cr_for_database(factor_common(1,1,1,1,1), size(factor_common,kind=4)) + if (size(scale_factor,kind=4) > 0) & + call bcast_all_cr_for_database(scale_factor(1,1,1,1), size(scale_factor,kind=4)) + call bcast_all_cr_for_database(factor_common_kappa(1,1,1,1,1), size(factor_common_kappa,kind=4)) + call bcast_all_cr_for_database(scale_factor_kappa(1,1,1,1), size(scale_factor_kappa,kind=4)) ! gets stress relaxation times tau_sigma, i.e. ! precalculates tau_sigma depending on period band (constant for all Q_mu), and diff --git a/src/specfem3D/prepare_gpu.f90 b/src/specfem3D/prepare_gpu.f90 index 6a5b8b3a3..b7136f58a 100644 --- a/src/specfem3D/prepare_gpu.f90 +++ b/src/specfem3D/prepare_gpu.f90 @@ -145,7 +145,7 @@ subroutine prepare_GPU() epsilondev_xx,epsilondev_yy,epsilondev_xy, & epsilondev_xz,epsilondev_yz, & ATTENUATION, & - size(R_xx), & + size(R_xx,kind=4), & R_xx,R_yy,R_xy,R_xz,R_yz, & factor_common, & R_trace,epsilondev_trace, & @@ -171,7 +171,7 @@ subroutine prepare_GPU() b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, & b_epsilondev_xz,b_epsilondev_yz, & b_epsilon_trace_over_3, & - size(R_xx), & + size(R_xx,kind=4), & b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, & b_R_trace,b_epsilondev_trace, & b_alphaval,b_betaval,b_gammaval, & diff --git a/src/specfem3D/read_forward_arrays.f90 b/src/specfem3D/read_forward_arrays.f90 index e76d90ab3..ca45d5033 100644 --- a/src/specfem3D/read_forward_arrays.f90 +++ b/src/specfem3D/read_forward_arrays.f90 @@ -110,11 +110,11 @@ subroutine read_forward_arrays() if (ATTENUATION) then call transfer_b_rmemory_to_device(Mesh_pointer, & b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, & - b_R_trace,size(b_R_xx)) + b_R_trace,size(b_R_xx,kind=4)) call transfer_b_strain_to_device(Mesh_pointer, & b_epsilondev_xx,b_epsilondev_yy,b_epsilondev_xy, & b_epsilondev_xz,b_epsilondev_yz, & - b_epsilondev_trace,size(b_epsilondev_xx)) + b_epsilondev_trace,size(b_epsilondev_xx,kind=4)) endif endif @@ -217,7 +217,7 @@ subroutine read_forward_arrays_undoatt() if (ATTENUATION) then call transfer_b_rmemory_to_device(Mesh_pointer, & b_R_xx,b_R_yy,b_R_xy,b_R_xz,b_R_yz, & - b_R_trace,size(b_R_xx)) + b_R_trace,size(b_R_xx,kind=4)) endif endif diff --git a/src/specfem3D/read_mesh_databases.F90 b/src/specfem3D/read_mesh_databases.F90 index be634c884..8f75106fe 100644 --- a/src/specfem3D/read_mesh_databases.F90 +++ b/src/specfem3D/read_mesh_databases.F90 @@ -192,35 +192,38 @@ subroutine read_mesh_databases() if (itest /= 9999) stop 'Error database read at position 1' endif +! note: the size(..) function returns either integer(kind=4) or integer(kind=8) +! depending on compiler flags (-mcmedium), thus adding a kind argument to have integer(kind=4) output + call bcast_all_i_for_database(NSPEC_AB, 1) call bcast_all_i_for_database(NGLOB_AB, 1) call bcast_all_i_for_database(NSPEC_IRREGULAR, 1) - call bcast_all_i_for_database(ibool(1,1,1,1), size(ibool)) - call bcast_all_i_for_database(irregular_element_number(1), size(irregular_element_number)) + call bcast_all_i_for_database(ibool(1,1,1,1), size(ibool,kind=4)) + call bcast_all_i_for_database(irregular_element_number(1), size(irregular_element_number,kind=4)) call bcast_all_cr_for_database(xix_regular, 1) call bcast_all_cr_for_database(jacobian_regular, 1) - call bcast_all_cr_for_database(xstore(1), size(xstore)) - call bcast_all_cr_for_database(ystore(1), size(ystore)) - call bcast_all_cr_for_database(zstore(1), size(zstore)) + call bcast_all_cr_for_database(xstore(1), size(xstore,kind=4)) + call bcast_all_cr_for_database(ystore(1), size(ystore,kind=4)) + call bcast_all_cr_for_database(zstore(1), size(zstore,kind=4)) - call bcast_all_cr_for_database(xixstore(1,1,1,1), size(xixstore)) - call bcast_all_cr_for_database(xiystore(1,1,1,1), size(xiystore)) - call bcast_all_cr_for_database(xizstore(1,1,1,1), size(xizstore)) - call bcast_all_cr_for_database(etaxstore(1,1,1,1), size(etaxstore)) - call bcast_all_cr_for_database(etaystore(1,1,1,1), size(etaystore)) - call bcast_all_cr_for_database(etazstore(1,1,1,1), size(etazstore)) - call bcast_all_cr_for_database(gammaxstore(1,1,1,1), size(gammaxstore)) - call bcast_all_cr_for_database(gammaystore(1,1,1,1), size(gammaystore)) - call bcast_all_cr_for_database(gammazstore(1,1,1,1), size(gammazstore)) - call bcast_all_cr_for_database(jacobianstore(1,1,1,1), size(jacobianstore)) + call bcast_all_cr_for_database(xixstore(1,1,1,1), size(xixstore,kind=4)) + call bcast_all_cr_for_database(xiystore(1,1,1,1), size(xiystore,kind=4)) + call bcast_all_cr_for_database(xizstore(1,1,1,1), size(xizstore,kind=4)) + call bcast_all_cr_for_database(etaxstore(1,1,1,1), size(etaxstore,kind=4)) + call bcast_all_cr_for_database(etaystore(1,1,1,1), size(etaystore,kind=4)) + call bcast_all_cr_for_database(etazstore(1,1,1,1), size(etazstore,kind=4)) + call bcast_all_cr_for_database(gammaxstore(1,1,1,1), size(gammaxstore,kind=4)) + call bcast_all_cr_for_database(gammaystore(1,1,1,1), size(gammaystore,kind=4)) + call bcast_all_cr_for_database(gammazstore(1,1,1,1), size(gammazstore,kind=4)) + call bcast_all_cr_for_database(jacobianstore(1,1,1,1), size(jacobianstore,kind=4)) - call bcast_all_cr_for_database(kappastore(1,1,1,1), size(kappastore)) - call bcast_all_cr_for_database(mustore(1,1,1,1), size(mustore)) + call bcast_all_cr_for_database(kappastore(1,1,1,1), size(kappastore,kind=4)) + call bcast_all_cr_for_database(mustore(1,1,1,1), size(mustore,kind=4)) - call bcast_all_l_for_database(ispec_is_acoustic(1), size(ispec_is_acoustic)) - call bcast_all_l_for_database(ispec_is_elastic(1), size(ispec_is_elastic)) - call bcast_all_l_for_database(ispec_is_poroelastic(1), size(ispec_is_poroelastic)) + call bcast_all_l_for_database(ispec_is_acoustic(1), size(ispec_is_acoustic,kind=4)) + call bcast_all_l_for_database(ispec_is_elastic(1), size(ispec_is_elastic,kind=4)) + call bcast_all_l_for_database(ispec_is_poroelastic(1), size(ispec_is_poroelastic,kind=4)) ! all processes will have acoustic_simulation set if any flag is .true. call any_all_l( ANY(ispec_is_acoustic), ACOUSTIC_SIMULATION ) @@ -277,7 +280,7 @@ subroutine read_mesh_databases() rmass_acoustic(:) = 0.0_CUSTOM_REAL if (I_should_read_the_database) read(IIN) rmass_acoustic - call bcast_all_cr_for_database(rmass_acoustic(1), size(rmass_acoustic)) + call bcast_all_cr_for_database(rmass_acoustic(1), size(rmass_acoustic,kind=4)) ! initializes mass matrix contribution allocate(rmassz_acoustic(NGLOB_AB),stat=ier) @@ -289,7 +292,7 @@ subroutine read_mesh_databases() ! rho array is needed for acoustic simulations but also for elastic simulations with CPML, ! read it in all cases (whether the simulation is acoustic, elastic, or acoustic/elastic) if (I_should_read_the_database) read(IIN) rhostore - call bcast_all_cr_for_database(rhostore(1,1,1,1), size(rhostore)) + call bcast_all_cr_for_database(rhostore(1,1,1,1), size(rhostore,kind=4)) ! elastic simulation if (ELASTIC_SIMULATION) then @@ -453,7 +456,7 @@ subroutine read_mesh_databases() read(IIN,iostat=ier) rmass if (ier /= 0) stop 'Error reading in array rmass' endif - call bcast_all_cr_for_database(rmass(1), size(rmass)) + call bcast_all_cr_for_database(rmass(1), size(rmass,kind=4)) if (APPROXIMATE_OCEAN_LOAD) then ! ocean mass matrix @@ -463,7 +466,7 @@ subroutine read_mesh_databases() rmass_ocean_load(:) = 0.0_CUSTOM_REAL if (I_should_read_the_database) read(IIN) rmass_ocean_load - call bcast_all_cr_for_database(rmass_ocean_load(1), size(rmass_ocean_load)) + call bcast_all_cr_for_database(rmass_ocean_load(1), size(rmass_ocean_load,kind=4)) else ! dummy allocation allocate(rmass_ocean_load(1),stat=ier) @@ -476,13 +479,13 @@ subroutine read_mesh_databases() read(IIN,iostat=ier) rho_vp if (ier /= 0) stop 'Error reading in array rho_vp' endif - call bcast_all_cr_for_database(rho_vp(1,1,1,1), size(rho_vp)) + call bcast_all_cr_for_database(rho_vp(1,1,1,1), size(rho_vp,kind=4)) if (I_should_read_the_database) then read(IIN,iostat=ier) rho_vs if (ier /= 0) stop 'Error reading in array rho_vs' endif - call bcast_all_cr_for_database(rho_vs(1,1,1,1), size(rho_vs)) + call bcast_all_cr_for_database(rho_vs(1,1,1,1), size(rho_vs,kind=4)) else ! no elastic attenuation & anisotropy @@ -604,17 +607,17 @@ subroutine read_mesh_databases() read(IIN) rho_vpII read(IIN) rho_vsI endif - call bcast_all_cr_for_database(rmass_solid_poroelastic(1), size(rmass_solid_poroelastic)) - call bcast_all_cr_for_database(rmass_fluid_poroelastic(1), size(rmass_fluid_poroelastic)) - call bcast_all_cr_for_database(rhoarraystore(1,1,1,1,1), size(rhoarraystore)) - call bcast_all_cr_for_database(kappaarraystore(1,1,1,1,1), size(kappaarraystore)) - call bcast_all_cr_for_database(etastore(1,1,1,1), size(etastore)) - call bcast_all_cr_for_database(tortstore(1,1,1,1), size(tortstore)) - call bcast_all_cr_for_database(permstore(1,1,1,1,1), size(permstore)) - call bcast_all_cr_for_database(phistore(1,1,1,1), size(phistore)) - call bcast_all_cr_for_database(rho_vpI(1,1,1,1), size(rho_vpI)) - call bcast_all_cr_for_database(rho_vpII(1,1,1,1), size(rho_vpII)) - call bcast_all_cr_for_database(rho_vsI(1,1,1,1), size(rho_vsI)) + call bcast_all_cr_for_database(rmass_solid_poroelastic(1), size(rmass_solid_poroelastic,kind=4)) + call bcast_all_cr_for_database(rmass_fluid_poroelastic(1), size(rmass_fluid_poroelastic,kind=4)) + call bcast_all_cr_for_database(rhoarraystore(1,1,1,1,1), size(rhoarraystore,kind=4)) + call bcast_all_cr_for_database(kappaarraystore(1,1,1,1,1), size(kappaarraystore,kind=4)) + call bcast_all_cr_for_database(etastore(1,1,1,1), size(etastore,kind=4)) + call bcast_all_cr_for_database(tortstore(1,1,1,1), size(tortstore,kind=4)) + call bcast_all_cr_for_database(permstore(1,1,1,1,1), size(permstore,kind=4)) + call bcast_all_cr_for_database(phistore(1,1,1,1), size(phistore,kind=4)) + call bcast_all_cr_for_database(rho_vpI(1,1,1,1), size(rho_vpI,kind=4)) + call bcast_all_cr_for_database(rho_vpII(1,1,1,1), size(rho_vpII,kind=4)) + call bcast_all_cr_for_database(rho_vsI(1,1,1,1), size(rho_vsI,kind=4)) else ! dummy allocations (needed for subroutine arguments) allocate(rhoarraystore(2,1,1,1,1), & @@ -716,18 +719,18 @@ subroutine read_mesh_databases() read(IIN) alpha_store_y read(IIN) alpha_store_z endif - call bcast_all_i_for_database(CPML_regions(1), size(CPML_regions)) - call bcast_all_i_for_database(CPML_to_spec(1), size(CPML_to_spec)) - call bcast_all_l_for_database(is_CPML(1), size(is_CPML)) - call bcast_all_cr_for_database(d_store_x(1,1,1,1), size(d_store_x)) - call bcast_all_cr_for_database(d_store_y(1,1,1,1), size(d_store_y)) - call bcast_all_cr_for_database(d_store_z(1,1,1,1), size(d_store_z)) - call bcast_all_cr_for_database(k_store_x(1,1,1,1), size(k_store_x)) - call bcast_all_cr_for_database(k_store_y(1,1,1,1), size(k_store_y)) - call bcast_all_cr_for_database(k_store_z(1,1,1,1), size(k_store_z)) - call bcast_all_cr_for_database(alpha_store_x(1,1,1,1), size(alpha_store_x)) - call bcast_all_cr_for_database(alpha_store_y(1,1,1,1), size(alpha_store_y)) - call bcast_all_cr_for_database(alpha_store_z(1,1,1,1), size(alpha_store_z)) + call bcast_all_i_for_database(CPML_regions(1), size(CPML_regions,kind=4)) + call bcast_all_i_for_database(CPML_to_spec(1), size(CPML_to_spec,kind=4)) + call bcast_all_l_for_database(is_CPML(1), size(is_CPML,kind=4)) + call bcast_all_cr_for_database(d_store_x(1,1,1,1), size(d_store_x,kind=4)) + call bcast_all_cr_for_database(d_store_y(1,1,1,1), size(d_store_y,kind=4)) + call bcast_all_cr_for_database(d_store_z(1,1,1,1), size(d_store_z,kind=4)) + call bcast_all_cr_for_database(k_store_x(1,1,1,1), size(k_store_x,kind=4)) + call bcast_all_cr_for_database(k_store_y(1,1,1,1), size(k_store_y,kind=4)) + call bcast_all_cr_for_database(k_store_z(1,1,1,1), size(k_store_z,kind=4)) + call bcast_all_cr_for_database(alpha_store_x(1,1,1,1), size(alpha_store_x,kind=4)) + call bcast_all_cr_for_database(alpha_store_y(1,1,1,1), size(alpha_store_y,kind=4)) + call bcast_all_cr_for_database(alpha_store_z(1,1,1,1), size(alpha_store_z,kind=4)) if ((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then if (I_should_read_the_database) read(IIN) nglob_interface_PML_acoustic @@ -743,7 +746,7 @@ subroutine read_mesh_databases() points_interface_PML_acoustic(:) = 0 if (I_should_read_the_database) read(IIN) points_interface_PML_acoustic - call bcast_all_i_for_database(points_interface_PML_acoustic(1), size(points_interface_PML_acoustic)) + call bcast_all_i_for_database(points_interface_PML_acoustic(1), size(points_interface_PML_acoustic,kind=4)) endif if (nglob_interface_PML_elastic > 0) then allocate(points_interface_PML_elastic(nglob_interface_PML_elastic),stat=ier) @@ -752,7 +755,7 @@ subroutine read_mesh_databases() points_interface_PML_elastic(:) = 0 if (I_should_read_the_database) read(IIN) points_interface_PML_elastic - call bcast_all_i_for_database(points_interface_PML_elastic(1), size(points_interface_PML_elastic)) + call bcast_all_i_for_database(points_interface_PML_elastic(1), size(points_interface_PML_elastic,kind=4)) endif endif endif @@ -786,10 +789,10 @@ subroutine read_mesh_databases() read(IIN) abs_boundary_jacobian2Dw read(IIN) abs_boundary_normal endif - call bcast_all_i_for_database(abs_boundary_ispec(1), size(abs_boundary_ispec)) - call bcast_all_i_for_database(abs_boundary_ijk(1,1,1), size(abs_boundary_ijk)) - call bcast_all_cr_for_database(abs_boundary_jacobian2Dw(1,1), size(abs_boundary_jacobian2Dw)) - call bcast_all_cr_for_database(abs_boundary_normal(1,1,1), size(abs_boundary_normal)) + call bcast_all_i_for_database(abs_boundary_ispec(1), size(abs_boundary_ispec,kind=4)) + call bcast_all_i_for_database(abs_boundary_ijk(1,1,1), size(abs_boundary_ijk,kind=4)) + call bcast_all_cr_for_database(abs_boundary_jacobian2Dw(1,1), size(abs_boundary_jacobian2Dw,kind=4)) + call bcast_all_cr_for_database(abs_boundary_normal(1,1,1), size(abs_boundary_normal,kind=4)) if (STACEY_ABSORBING_CONDITIONS .and. (.not. PML_CONDITIONS)) then ! store mass matrix contributions @@ -799,13 +802,13 @@ subroutine read_mesh_databases() read(IIN) rmassy read(IIN) rmassz endif - call bcast_all_cr_for_database(rmassx(1), size(rmassx)) - call bcast_all_cr_for_database(rmassy(1), size(rmassy)) - call bcast_all_cr_for_database(rmassz(1), size(rmassz)) + call bcast_all_cr_for_database(rmassx(1), size(rmassx,kind=4)) + call bcast_all_cr_for_database(rmassy(1), size(rmassy,kind=4)) + call bcast_all_cr_for_database(rmassz(1), size(rmassz,kind=4)) endif if (ACOUSTIC_SIMULATION) then if (I_should_read_the_database) read(IIN) rmassz_acoustic - call bcast_all_cr_for_database(rmassz_acoustic(1), size(rmassz_acoustic)) + call bcast_all_cr_for_database(rmassz_acoustic(1), size(rmassz_acoustic,kind=4)) endif endif endif @@ -857,12 +860,12 @@ subroutine read_mesh_databases() if (nspec2D_bottom > 0) read(IIN) ibelm_bottom if (nspec2D_top > 0) read(IIN) ibelm_top endif - if (nspec2D_xmin > 0) call bcast_all_i_for_database(ibelm_xmin(1), size(ibelm_xmin)) - if (nspec2D_xmax > 0) call bcast_all_i_for_database(ibelm_xmax(1), size(ibelm_xmax)) - if (nspec2D_ymin > 0) call bcast_all_i_for_database(ibelm_ymin(1), size(ibelm_ymin)) - if (nspec2D_ymax > 0) call bcast_all_i_for_database(ibelm_ymax(1), size(ibelm_ymax)) - if (nspec2D_bottom > 0) call bcast_all_i_for_database(ibelm_bottom(1), size(ibelm_bottom)) - if (nspec2D_top > 0) call bcast_all_i_for_database(ibelm_top(1), size(ibelm_top)) + if (nspec2D_xmin > 0) call bcast_all_i_for_database(ibelm_xmin(1), size(ibelm_xmin,kind=4)) + if (nspec2D_xmax > 0) call bcast_all_i_for_database(ibelm_xmax(1), size(ibelm_xmax,kind=4)) + if (nspec2D_ymin > 0) call bcast_all_i_for_database(ibelm_ymin(1), size(ibelm_ymin,kind=4)) + if (nspec2D_ymax > 0) call bcast_all_i_for_database(ibelm_ymax(1), size(ibelm_ymax,kind=4)) + if (nspec2D_bottom > 0) call bcast_all_i_for_database(ibelm_bottom(1), size(ibelm_bottom,kind=4)) + if (nspec2D_top > 0) call bcast_all_i_for_database(ibelm_top(1), size(ibelm_top,kind=4)) ! free surface if (I_should_read_the_database) read(IIN) num_free_surface_faces @@ -887,10 +890,10 @@ subroutine read_mesh_databases() read(IIN) free_surface_jacobian2Dw read(IIN) free_surface_normal endif - call bcast_all_i_for_database(free_surface_ispec(1), size(free_surface_ispec)) - call bcast_all_i_for_database(free_surface_ijk(1,1,1), size(free_surface_ijk)) - call bcast_all_cr_for_database(free_surface_jacobian2Dw(1,1), size(free_surface_jacobian2Dw)) - call bcast_all_cr_for_database(free_surface_normal(1,1,1), size(free_surface_normal)) + call bcast_all_i_for_database(free_surface_ispec(1), size(free_surface_ispec,kind=4)) + call bcast_all_i_for_database(free_surface_ijk(1,1,1), size(free_surface_ijk,kind=4)) + call bcast_all_cr_for_database(free_surface_jacobian2Dw(1,1), size(free_surface_jacobian2Dw,kind=4)) + call bcast_all_cr_for_database(free_surface_normal(1,1,1), size(free_surface_normal,kind=4)) endif ! acoustic-elastic coupling surface @@ -916,10 +919,10 @@ subroutine read_mesh_databases() read(IIN) coupling_ac_el_jacobian2Dw read(IIN) coupling_ac_el_normal endif - call bcast_all_i_for_database(coupling_ac_el_ispec(1), size(coupling_ac_el_ispec)) - call bcast_all_i_for_database(coupling_ac_el_ijk(1,1,1), size(coupling_ac_el_ijk)) - call bcast_all_cr_for_database(coupling_ac_el_jacobian2Dw(1,1), size(coupling_ac_el_jacobian2Dw)) - call bcast_all_cr_for_database(coupling_ac_el_normal(1,1,1), size(coupling_ac_el_normal)) + call bcast_all_i_for_database(coupling_ac_el_ispec(1), size(coupling_ac_el_ispec,kind=4)) + call bcast_all_i_for_database(coupling_ac_el_ijk(1,1,1), size(coupling_ac_el_ijk,kind=4)) + call bcast_all_cr_for_database(coupling_ac_el_jacobian2Dw(1,1), size(coupling_ac_el_jacobian2Dw,kind=4)) + call bcast_all_cr_for_database(coupling_ac_el_normal(1,1,1), size(coupling_ac_el_normal,kind=4)) endif ! acoustic-poroelastic coupling surface @@ -945,10 +948,10 @@ subroutine read_mesh_databases() read(IIN) coupling_ac_po_jacobian2Dw read(IIN) coupling_ac_po_normal endif - call bcast_all_i_for_database(coupling_ac_po_ispec(1), size(coupling_ac_po_ispec)) - call bcast_all_i_for_database(coupling_ac_po_ijk(1,1,1), size(coupling_ac_po_ijk)) - call bcast_all_cr_for_database(coupling_ac_po_jacobian2Dw(1,1), size(coupling_ac_po_jacobian2Dw)) - call bcast_all_cr_for_database(coupling_ac_po_normal(1,1,1), size(coupling_ac_po_normal)) + call bcast_all_i_for_database(coupling_ac_po_ispec(1), size(coupling_ac_po_ispec,kind=4)) + call bcast_all_i_for_database(coupling_ac_po_ijk(1,1,1), size(coupling_ac_po_ijk,kind=4)) + call bcast_all_cr_for_database(coupling_ac_po_jacobian2Dw(1,1), size(coupling_ac_po_jacobian2Dw,kind=4)) + call bcast_all_cr_for_database(coupling_ac_po_normal(1,1,1), size(coupling_ac_po_normal,kind=4)) endif ! elastic-poroelastic coupling surface @@ -981,12 +984,12 @@ subroutine read_mesh_databases() read(IIN) coupling_el_po_jacobian2Dw read(IIN) coupling_el_po_normal endif - call bcast_all_i_for_database(coupling_el_po_ispec(1), size(coupling_el_po_ispec)) - call bcast_all_i_for_database(coupling_po_el_ispec(1), size(coupling_po_el_ispec)) - call bcast_all_i_for_database(coupling_el_po_ijk(1,1,1), size(coupling_el_po_ijk)) - call bcast_all_i_for_database(coupling_po_el_ijk(1,1,1), size(coupling_po_el_ijk)) - call bcast_all_cr_for_database(coupling_el_po_jacobian2Dw(1,1), size(coupling_el_po_jacobian2Dw)) - call bcast_all_cr_for_database(coupling_el_po_normal(1,1,1), size(coupling_el_po_normal)) + call bcast_all_i_for_database(coupling_el_po_ispec(1), size(coupling_el_po_ispec,kind=4)) + call bcast_all_i_for_database(coupling_po_el_ispec(1), size(coupling_po_el_ispec,kind=4)) + call bcast_all_i_for_database(coupling_el_po_ijk(1,1,1), size(coupling_el_po_ijk,kind=4)) + call bcast_all_i_for_database(coupling_po_el_ijk(1,1,1), size(coupling_po_el_ijk,kind=4)) + call bcast_all_cr_for_database(coupling_el_po_jacobian2Dw(1,1), size(coupling_el_po_jacobian2Dw,kind=4)) + call bcast_all_cr_for_database(coupling_el_po_normal(1,1,1), size(coupling_el_po_normal,kind=4)) endif ! checks i/o so far @@ -1021,9 +1024,9 @@ subroutine read_mesh_databases() read(IIN) nibool_interfaces_ext_mesh read(IIN) ibool_interfaces_ext_mesh endif - call bcast_all_i_for_database(my_neighbors_ext_mesh(1), size(my_neighbors_ext_mesh)) - call bcast_all_i_for_database(nibool_interfaces_ext_mesh(1), size(nibool_interfaces_ext_mesh)) - call bcast_all_i_for_database(ibool_interfaces_ext_mesh(1,1), size(ibool_interfaces_ext_mesh)) + call bcast_all_i_for_database(my_neighbors_ext_mesh(1), size(my_neighbors_ext_mesh,kind=4)) + call bcast_all_i_for_database(nibool_interfaces_ext_mesh(1), size(nibool_interfaces_ext_mesh,kind=4)) + call bcast_all_i_for_database(ibool_interfaces_ext_mesh(1,1), size(ibool_interfaces_ext_mesh,kind=4)) else ! no interfaces max_nibool_interfaces_ext_mesh = 0 @@ -1067,27 +1070,27 @@ subroutine read_mesh_databases() read(IIN) c56store read(IIN) c66store endif - call bcast_all_cr_for_database(c11store(1,1,1,1), size(c11store)) - call bcast_all_cr_for_database(c12store(1,1,1,1), size(c12store)) - call bcast_all_cr_for_database(c13store(1,1,1,1), size(c13store)) - call bcast_all_cr_for_database(c14store(1,1,1,1), size(c14store)) - call bcast_all_cr_for_database(c15store(1,1,1,1), size(c15store)) - call bcast_all_cr_for_database(c16store(1,1,1,1), size(c16store)) - call bcast_all_cr_for_database(c22store(1,1,1,1), size(c22store)) - call bcast_all_cr_for_database(c23store(1,1,1,1), size(c23store)) - call bcast_all_cr_for_database(c24store(1,1,1,1), size(c24store)) - call bcast_all_cr_for_database(c25store(1,1,1,1), size(c25store)) - call bcast_all_cr_for_database(c26store(1,1,1,1), size(c26store)) - call bcast_all_cr_for_database(c33store(1,1,1,1), size(c33store)) - call bcast_all_cr_for_database(c34store(1,1,1,1), size(c34store)) - call bcast_all_cr_for_database(c35store(1,1,1,1), size(c35store)) - call bcast_all_cr_for_database(c36store(1,1,1,1), size(c36store)) - call bcast_all_cr_for_database(c44store(1,1,1,1), size(c44store)) - call bcast_all_cr_for_database(c45store(1,1,1,1), size(c45store)) - call bcast_all_cr_for_database(c46store(1,1,1,1), size(c46store)) - call bcast_all_cr_for_database(c55store(1,1,1,1), size(c55store)) - call bcast_all_cr_for_database(c56store(1,1,1,1), size(c56store)) - call bcast_all_cr_for_database(c66store(1,1,1,1), size(c66store)) + call bcast_all_cr_for_database(c11store(1,1,1,1), size(c11store,kind=4)) + call bcast_all_cr_for_database(c12store(1,1,1,1), size(c12store,kind=4)) + call bcast_all_cr_for_database(c13store(1,1,1,1), size(c13store,kind=4)) + call bcast_all_cr_for_database(c14store(1,1,1,1), size(c14store,kind=4)) + call bcast_all_cr_for_database(c15store(1,1,1,1), size(c15store,kind=4)) + call bcast_all_cr_for_database(c16store(1,1,1,1), size(c16store,kind=4)) + call bcast_all_cr_for_database(c22store(1,1,1,1), size(c22store,kind=4)) + call bcast_all_cr_for_database(c23store(1,1,1,1), size(c23store,kind=4)) + call bcast_all_cr_for_database(c24store(1,1,1,1), size(c24store,kind=4)) + call bcast_all_cr_for_database(c25store(1,1,1,1), size(c25store,kind=4)) + call bcast_all_cr_for_database(c26store(1,1,1,1), size(c26store,kind=4)) + call bcast_all_cr_for_database(c33store(1,1,1,1), size(c33store,kind=4)) + call bcast_all_cr_for_database(c34store(1,1,1,1), size(c34store,kind=4)) + call bcast_all_cr_for_database(c35store(1,1,1,1), size(c35store,kind=4)) + call bcast_all_cr_for_database(c36store(1,1,1,1), size(c36store,kind=4)) + call bcast_all_cr_for_database(c44store(1,1,1,1), size(c44store,kind=4)) + call bcast_all_cr_for_database(c45store(1,1,1,1), size(c45store,kind=4)) + call bcast_all_cr_for_database(c46store(1,1,1,1), size(c46store,kind=4)) + call bcast_all_cr_for_database(c55store(1,1,1,1), size(c55store,kind=4)) + call bcast_all_cr_for_database(c56store(1,1,1,1), size(c56store,kind=4)) + call bcast_all_cr_for_database(c66store(1,1,1,1), size(c66store,kind=4)) endif ! inner / outer elements @@ -1097,7 +1100,7 @@ subroutine read_mesh_databases() ispec_is_inner(:) = .false. if (I_should_read_the_database) read(IIN) ispec_is_inner - call bcast_all_l_for_database(ispec_is_inner(1), size(ispec_is_inner)) + call bcast_all_l_for_database(ispec_is_inner(1), size(ispec_is_inner,kind=4)) if (ACOUSTIC_SIMULATION) then if (I_should_read_the_database) then @@ -1116,7 +1119,7 @@ subroutine read_mesh_databases() if (num_phase_ispec_acoustic > 0) then if (I_should_read_the_database) read(IIN) phase_ispec_inner_acoustic - call bcast_all_i_for_database(phase_ispec_inner_acoustic(1,1), size(phase_ispec_inner_acoustic)) + call bcast_all_i_for_database(phase_ispec_inner_acoustic(1,1), size(phase_ispec_inner_acoustic,kind=4)) endif endif @@ -1137,7 +1140,7 @@ subroutine read_mesh_databases() if (num_phase_ispec_elastic > 0) then if (I_should_read_the_database) read(IIN) phase_ispec_inner_elastic - call bcast_all_i_for_database(phase_ispec_inner_elastic(1,1), size(phase_ispec_inner_elastic)) + call bcast_all_i_for_database(phase_ispec_inner_elastic(1,1), size(phase_ispec_inner_elastic,kind=4)) endif endif @@ -1158,7 +1161,7 @@ subroutine read_mesh_databases() if (num_phase_ispec_poroelastic > 0) then if (I_should_read_the_database) read(IIN) phase_ispec_inner_poroelastic - call bcast_all_i_for_database(phase_ispec_inner_poroelastic(1,1), size(phase_ispec_inner_poroelastic)) + call bcast_all_i_for_database(phase_ispec_inner_poroelastic(1,1), size(phase_ispec_inner_poroelastic,kind=4)) endif endif @@ -1176,7 +1179,7 @@ subroutine read_mesh_databases() num_elem_colors_acoustic(:) = 0 if (I_should_read_the_database) read(IIN) num_elem_colors_acoustic - call bcast_all_i_for_database(num_elem_colors_acoustic(1), size(num_elem_colors_acoustic)) + call bcast_all_i_for_database(num_elem_colors_acoustic(1), size(num_elem_colors_acoustic,kind=4)) endif ! elastic domain colors if (ELASTIC_SIMULATION) then @@ -1190,7 +1193,7 @@ subroutine read_mesh_databases() num_elem_colors_elastic(:) = 0 if (I_should_read_the_database) read(IIN) num_elem_colors_elastic - call bcast_all_i_for_database(num_elem_colors_elastic(1), size(num_elem_colors_elastic)) + call bcast_all_i_for_database(num_elem_colors_elastic(1), size(num_elem_colors_elastic,kind=4)) endif else ! allocates dummy arrays @@ -1228,8 +1231,8 @@ subroutine read_mesh_databases() read(IIN) iglob_is_surface_external_mesh endif call bcast_all_i_for_database(nfaces_surface, 1) - call bcast_all_l_for_database(ispec_is_surface_external_mesh(1), size(ispec_is_surface_external_mesh)) - call bcast_all_l_for_database(iglob_is_surface_external_mesh(1), size(iglob_is_surface_external_mesh)) + call bcast_all_l_for_database(ispec_is_surface_external_mesh(1), size(ispec_is_surface_external_mesh,kind=4)) + call bcast_all_l_for_database(iglob_is_surface_external_mesh(1), size(iglob_is_surface_external_mesh,kind=4)) ! for mesh adjacency if (I_should_read_the_database) then @@ -1247,8 +1250,8 @@ subroutine read_mesh_databases() read(IIN) neighbors_xadj read(IIN) neighbors_adjncy endif - call bcast_all_i_for_database(neighbors_xadj(1), size(neighbors_xadj)) - call bcast_all_i_for_database(neighbors_adjncy(1), size(neighbors_adjncy)) + call bcast_all_i_for_database(neighbors_xadj(1), size(neighbors_xadj,kind=4)) + call bcast_all_i_for_database(neighbors_adjncy(1), size(neighbors_adjncy,kind=4)) ! checks i/o so far if (I_should_read_the_database) then @@ -1484,10 +1487,10 @@ subroutine read_mesh_databases_moho() read(IIN) ijk_moho_top read(IIN) ijk_moho_bot endif - call bcast_all_i_for_database(ibelm_moho_top(1), size(ibelm_moho_top)) - call bcast_all_i_for_database(ibelm_moho_bot(1), size(ibelm_moho_bot)) - call bcast_all_i_for_database(ijk_moho_top(1,1,1), size(ijk_moho_top)) - call bcast_all_i_for_database(ijk_moho_bot(1,1,1), size(ijk_moho_bot)) + call bcast_all_i_for_database(ibelm_moho_top(1), size(ibelm_moho_top,kind=4)) + call bcast_all_i_for_database(ibelm_moho_bot(1), size(ibelm_moho_bot,kind=4)) + call bcast_all_i_for_database(ijk_moho_top(1,1,1), size(ijk_moho_top,kind=4)) + call bcast_all_i_for_database(ijk_moho_bot(1,1,1), size(ijk_moho_bot,kind=4)) if (I_should_read_the_database) close(IIN) @@ -1505,8 +1508,8 @@ subroutine read_mesh_databases_moho() read(IIN) normal_moho_top read(IIN) normal_moho_bot endif - call bcast_all_cr_for_database(normal_moho_top(1,1,1), size(normal_moho_top)) - call bcast_all_cr_for_database(normal_moho_bot(1,1,1), size(normal_moho_bot)) + call bcast_all_cr_for_database(normal_moho_top(1,1,1), size(normal_moho_top,kind=4)) + call bcast_all_cr_for_database(normal_moho_bot(1,1,1), size(normal_moho_bot,kind=4)) if (I_should_read_the_database) close(IIN) @@ -1524,8 +1527,8 @@ subroutine read_mesh_databases_moho() read(IIN) is_moho_top read(IIN) is_moho_bot endif - call bcast_all_l_for_database(is_moho_top(1), size(is_moho_top)) - call bcast_all_l_for_database(is_moho_bot(1), size(is_moho_bot)) + call bcast_all_l_for_database(is_moho_top(1), size(is_moho_top,kind=4)) + call bcast_all_l_for_database(is_moho_bot(1), size(is_moho_bot,kind=4)) if (I_should_read_the_database) close(IIN) diff --git a/src/specfem3D/read_mesh_databases_hdf5.F90 b/src/specfem3D/read_mesh_databases_hdf5.F90 index 0b82b5556..817258550 100644 --- a/src/specfem3D/read_mesh_databases_hdf5.F90 +++ b/src/specfem3D/read_mesh_databases_hdf5.F90 @@ -267,38 +267,41 @@ subroutine read_mesh_databases_hdf5() call h5_read_dataset_collect_hyperslab(dsetname,mustore,(/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL) endif +! note: the size(..) function returns either integer(kind=4) or integer(kind=8) +! depending on compiler flags (-mcmedium), thus adding a kind argument to have integer(kind=4) output + call bcast_all_i_for_database(NSPEC_AB, 1) call bcast_all_i_for_database(NGLOB_AB, 1) call bcast_all_i_for_database(NSPEC_IRREGULAR, 1) - call bcast_all_i_for_database(ibool(1,1,1,1), size(ibool)) - call bcast_all_i_for_database(irregular_element_number(1), size(irregular_element_number)) + call bcast_all_i_for_database(ibool(1,1,1,1), size(ibool,kind=4)) + call bcast_all_i_for_database(irregular_element_number(1), size(irregular_element_number,kind=4)) call bcast_all_cr_for_database(xix_regular, 1) call bcast_all_cr_for_database(jacobian_regular, 1) - if (size(xstore) > 0) call bcast_all_cr_for_database(xstore(1), size(xstore)) - if (size(ystore) > 0) call bcast_all_cr_for_database(ystore(1), size(ystore)) - if (size(zstore) > 0) call bcast_all_cr_for_database(zstore(1), size(zstore)) + if (size(xstore) > 0) call bcast_all_cr_for_database(xstore(1), size(xstore,kind=4)) + if (size(ystore) > 0) call bcast_all_cr_for_database(ystore(1), size(ystore,kind=4)) + if (size(zstore) > 0) call bcast_all_cr_for_database(zstore(1), size(zstore,kind=4)) - call bcast_all_cr_for_database(xixstore(1,1,1,1), size(xixstore)) - call bcast_all_cr_for_database(xiystore(1,1,1,1), size(xiystore)) - call bcast_all_cr_for_database(xizstore(1,1,1,1), size(xizstore)) - call bcast_all_cr_for_database(etaxstore(1,1,1,1), size(etaxstore)) - call bcast_all_cr_for_database(etaystore(1,1,1,1), size(etaystore)) - call bcast_all_cr_for_database(etazstore(1,1,1,1), size(etazstore)) - call bcast_all_cr_for_database(gammaxstore(1,1,1,1), size(gammaxstore)) - call bcast_all_cr_for_database(gammaystore(1,1,1,1), size(gammaystore)) - call bcast_all_cr_for_database(gammazstore(1,1,1,1), size(gammazstore)) - call bcast_all_cr_for_database(jacobianstore(1,1,1,1), size(jacobianstore)) + call bcast_all_cr_for_database(xixstore(1,1,1,1), size(xixstore,kind=4)) + call bcast_all_cr_for_database(xiystore(1,1,1,1), size(xiystore,kind=4)) + call bcast_all_cr_for_database(xizstore(1,1,1,1), size(xizstore,kind=4)) + call bcast_all_cr_for_database(etaxstore(1,1,1,1), size(etaxstore,kind=4)) + call bcast_all_cr_for_database(etaystore(1,1,1,1), size(etaystore,kind=4)) + call bcast_all_cr_for_database(etazstore(1,1,1,1), size(etazstore,kind=4)) + call bcast_all_cr_for_database(gammaxstore(1,1,1,1), size(gammaxstore,kind=4)) + call bcast_all_cr_for_database(gammaystore(1,1,1,1), size(gammaystore,kind=4)) + call bcast_all_cr_for_database(gammazstore(1,1,1,1), size(gammazstore,kind=4)) + call bcast_all_cr_for_database(jacobianstore(1,1,1,1), size(jacobianstore,kind=4)) - call bcast_all_cr_for_database(kappastore(1,1,1,1), size(kappastore)) - call bcast_all_cr_for_database(mustore(1,1,1,1), size(mustore)) + call bcast_all_cr_for_database(kappastore(1,1,1,1), size(kappastore,kind=4)) + call bcast_all_cr_for_database(mustore(1,1,1,1), size(mustore,kind=4)) if (size(ispec_is_acoustic) > 0) & - call bcast_all_l_for_database(ispec_is_acoustic(1), size(ispec_is_acoustic)) + call bcast_all_l_for_database(ispec_is_acoustic(1), size(ispec_is_acoustic,kind=4)) if (size(ispec_is_elastic) > 0) & - call bcast_all_l_for_database(ispec_is_elastic(1), size(ispec_is_elastic)) + call bcast_all_l_for_database(ispec_is_elastic(1), size(ispec_is_elastic,kind=4)) if (size(ispec_is_poroelastic) > 0) & - call bcast_all_l_for_database(ispec_is_poroelastic(1), size(ispec_is_poroelastic)) + call bcast_all_l_for_database(ispec_is_poroelastic(1), size(ispec_is_poroelastic,kind=4)) ! all processes will have acoustic_simulation set if any flag is .true. call any_all_l( ANY(ispec_is_acoustic), ACOUSTIC_SIMULATION ) @@ -357,7 +360,7 @@ subroutine read_mesh_databases_hdf5() dsetname = "rmass_acoustic" call h5_read_dataset_collect_hyperslab(dsetname, rmass_acoustic,(/sum(offset_nglob(0:myrank-1))/),H5_COL) endif - if (size(rmass_acoustic) > 0) call bcast_all_cr_for_database(rmass_acoustic(1), size(rmass_acoustic)) + if (size(rmass_acoustic) > 0) call bcast_all_cr_for_database(rmass_acoustic(1), size(rmass_acoustic,kind=4)) ! initializes mass matrix contribution allocate(rmassz_acoustic(NGLOB_AB),stat=ier) @@ -373,7 +376,7 @@ subroutine read_mesh_databases_hdf5() !call h5_read_dataset_p(dsetname, rhostore) call h5_read_dataset_collect_hyperslab(dsetname, rhostore, (/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL) endif - call bcast_all_cr_for_database(rhostore(1,1,1,1), size(rhostore)) + call bcast_all_cr_for_database(rhostore(1,1,1,1), size(rhostore,kind=4)) ! elastic simulation if (ELASTIC_SIMULATION) then @@ -537,7 +540,7 @@ subroutine read_mesh_databases_hdf5() dsetname = "rmass" call h5_read_dataset_collect_hyperslab(dsetname, rmass, (/sum(offset_nglob(0:myrank-1))/), H5_COL) endif - call bcast_all_cr_for_database(rmass(1), size(rmass)) + call bcast_all_cr_for_database(rmass(1), size(rmass,kind=4)) if (ier /= 0) stop 'Error reading in array rmass' if (APPROXIMATE_OCEAN_LOAD) then @@ -552,7 +555,7 @@ subroutine read_mesh_databases_hdf5() dsetname = "rmass_ocean_load" call h5_read_dataset_collect_hyperslab(dsetname, rmass_ocean_load, (/sum(offset_nglob_ocean(0:myrank-1))/),H5_COL) endif - if (size(rmass_ocean_load) > 0) call bcast_all_cr_for_database(rmass_ocean_load(1), size(rmass_ocean_load)) + if (size(rmass_ocean_load) > 0) call bcast_all_cr_for_database(rmass_ocean_load(1), size(rmass_ocean_load,kind=4)) else ! dummy allocation allocate(rmass_ocean_load(1),stat=ier) @@ -565,12 +568,12 @@ subroutine read_mesh_databases_hdf5() dsetname = "rho_vp" call h5_read_dataset_collect_hyperslab(dsetname, rho_vp, (/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL) endif - if (size(rho_vp) > 0) call bcast_all_cr_for_database(rho_vp(1,1,1,1), size(rho_vp)) + if (size(rho_vp) > 0) call bcast_all_cr_for_database(rho_vp(1,1,1,1), size(rho_vp,kind=4)) if (I_should_read_the_database) then dsetname = "rho_vs" call h5_read_dataset_collect_hyperslab(dsetname, rho_vs, (/0,0,0,sum(offset_nspec(0:myrank-1))/),H5_COL) endif - if (size(rho_vs) > 0) call bcast_all_cr_for_database(rho_vs(1,1,1,1), size(rho_vs)) + if (size(rho_vs) > 0) call bcast_all_cr_for_database(rho_vs(1,1,1,1), size(rho_vs,kind=4)) else ! no elastic attenuation & anisotropy @@ -706,27 +709,27 @@ subroutine read_mesh_databases_hdf5() call h5_read_dataset_collect_hyperslab(dsetname, rho_vsI, (/0,0,0,sum(offset_nspecporo(0:myrank-1))/),H5_COL) endif if (size(rmass_solid_poroelastic) > 0) & - call bcast_all_cr_for_database(rmass_solid_poroelastic(1), size(rmass_solid_poroelastic)) + call bcast_all_cr_for_database(rmass_solid_poroelastic(1), size(rmass_solid_poroelastic,kind=4)) if (size(rmass_fluid_poroelastic) > 0) & - call bcast_all_cr_for_database(rmass_fluid_poroelastic(1), size(rmass_fluid_poroelastic)) + call bcast_all_cr_for_database(rmass_fluid_poroelastic(1), size(rmass_fluid_poroelastic,kind=4)) if (size(rhoarraystore) > 0) & - call bcast_all_cr_for_database(rhoarraystore(1,1,1,1,1), size(rhoarraystore)) + call bcast_all_cr_for_database(rhoarraystore(1,1,1,1,1), size(rhoarraystore,kind=4)) if (size(kappaarraystore) > 0) & - call bcast_all_cr_for_database(kappaarraystore(1,1,1,1,1), size(kappaarraystore)) + call bcast_all_cr_for_database(kappaarraystore(1,1,1,1,1), size(kappaarraystore,kind=4)) if (size(etastore) > 0) & - call bcast_all_cr_for_database(etastore(1,1,1,1), size(etastore)) + call bcast_all_cr_for_database(etastore(1,1,1,1), size(etastore,kind=4)) if (size(tortstore) > 0) & - call bcast_all_cr_for_database(tortstore(1,1,1,1), size(tortstore)) + call bcast_all_cr_for_database(tortstore(1,1,1,1), size(tortstore,kind=4)) if (size(permstore) > 0) & - call bcast_all_cr_for_database(permstore(1,1,1,1,1), size(permstore)) + call bcast_all_cr_for_database(permstore(1,1,1,1,1), size(permstore,kind=4)) if (size(phistore) > 0) & - call bcast_all_cr_for_database(phistore(1,1,1,1), size(phistore)) + call bcast_all_cr_for_database(phistore(1,1,1,1), size(phistore,kind=4)) if (size(rho_vpI) > 0) & - call bcast_all_cr_for_database(rho_vpI(1,1,1,1), size(rho_vpI)) + call bcast_all_cr_for_database(rho_vpI(1,1,1,1), size(rho_vpI,kind=4)) if (size(rho_vpII) > 0) & - call bcast_all_cr_for_database(rho_vpII(1,1,1,1), size(rho_vpII)) + call bcast_all_cr_for_database(rho_vpII(1,1,1,1), size(rho_vpII,kind=4)) if (size(rho_vsI) > 0) & - call bcast_all_cr_for_database(rho_vsI(1,1,1,1), size(rho_vsI)) + call bcast_all_cr_for_database(rho_vsI(1,1,1,1), size(rho_vsI,kind=4)) else ! dummy allocations (needed for subroutine arguments) allocate(rhoarraystore(2,1,1,1,1), & @@ -846,18 +849,18 @@ subroutine read_mesh_databases_hdf5() call h5_read_dataset_collect_hyperslab(dsetname, alpha_store_z, (/0,0,0,sum(offset_nspeccpml(0:myrank-1))/),H5_COL) endif - if (size(CPML_regions) > 0) call bcast_all_i_for_database(CPML_regions(1), size(CPML_regions)) - if (size(CPML_to_spec) > 0) call bcast_all_i_for_database(CPML_to_spec(1), size(CPML_to_spec)) - if (size(is_CPML) > 0) call bcast_all_l_for_database(is_CPML(1), size(is_CPML)) - if (size(d_store_x) > 0) call bcast_all_cr_for_database(d_store_x(1,1,1,1), size(d_store_x)) - if (size(d_store_y) > 0) call bcast_all_cr_for_database(d_store_y(1,1,1,1), size(d_store_y)) - if (size(d_store_z) > 0) call bcast_all_cr_for_database(d_store_z(1,1,1,1), size(d_store_z)) - if (size(k_store_x) > 0) call bcast_all_cr_for_database(k_store_x(1,1,1,1), size(k_store_x)) - if (size(k_store_y) > 0) call bcast_all_cr_for_database(k_store_y(1,1,1,1), size(k_store_y)) - if (size(k_store_z) > 0) call bcast_all_cr_for_database(k_store_z(1,1,1,1), size(k_store_z)) - if (size(alpha_store_x) > 0) call bcast_all_cr_for_database(alpha_store_x(1,1,1,1), size(alpha_store_x)) - if (size(alpha_store_y) > 0) call bcast_all_cr_for_database(alpha_store_y(1,1,1,1), size(alpha_store_y)) - if (size(alpha_store_z) > 0) call bcast_all_cr_for_database(alpha_store_z(1,1,1,1), size(alpha_store_z)) + if (size(CPML_regions) > 0) call bcast_all_i_for_database(CPML_regions(1), size(CPML_regions,kind=4)) + if (size(CPML_to_spec) > 0) call bcast_all_i_for_database(CPML_to_spec(1), size(CPML_to_spec,kind=4)) + if (size(is_CPML) > 0) call bcast_all_l_for_database(is_CPML(1), size(is_CPML,kind=4)) + if (size(d_store_x) > 0) call bcast_all_cr_for_database(d_store_x(1,1,1,1), size(d_store_x,kind=4)) + if (size(d_store_y) > 0) call bcast_all_cr_for_database(d_store_y(1,1,1,1), size(d_store_y,kind=4)) + if (size(d_store_z) > 0) call bcast_all_cr_for_database(d_store_z(1,1,1,1), size(d_store_z,kind=4)) + if (size(k_store_x) > 0) call bcast_all_cr_for_database(k_store_x(1,1,1,1), size(k_store_x,kind=4)) + if (size(k_store_y) > 0) call bcast_all_cr_for_database(k_store_y(1,1,1,1), size(k_store_y,kind=4)) + if (size(k_store_z) > 0) call bcast_all_cr_for_database(k_store_z(1,1,1,1), size(k_store_z,kind=4)) + if (size(alpha_store_x) > 0) call bcast_all_cr_for_database(alpha_store_x(1,1,1,1), size(alpha_store_x,kind=4)) + if (size(alpha_store_y) > 0) call bcast_all_cr_for_database(alpha_store_y(1,1,1,1), size(alpha_store_y,kind=4)) + if (size(alpha_store_z) > 0) call bcast_all_cr_for_database(alpha_store_z(1,1,1,1), size(alpha_store_z,kind=4)) if ((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then ! acoustic @@ -889,7 +892,7 @@ subroutine read_mesh_databases_hdf5() points_interface_PML_acoustic, (/sum(offset_nglob_interface_PML_acoustic(0:myrank-1))/),H5_COL) endif if (size(points_interface_PML_acoustic) > 0) & - call bcast_all_i_for_database(points_interface_PML_acoustic(1), size(points_interface_PML_acoustic)) + call bcast_all_i_for_database(points_interface_PML_acoustic(1), size(points_interface_PML_acoustic,kind=4)) endif ! elastic allocation if (sum(offset_nglob_interface_PML_elastic) > 0) then @@ -904,7 +907,7 @@ subroutine read_mesh_databases_hdf5() points_interface_PML_acoustic, (/sum(offset_nglob_interface_PML_elastic(0:myrank-1))/),H5_COL) endif if (size(points_interface_PML_elastic) > 0) & - call bcast_all_i_for_database(points_interface_PML_elastic(1), size(points_interface_PML_elastic)) + call bcast_all_i_for_database(points_interface_PML_elastic(1), size(points_interface_PML_elastic,kind=4)) endif endif endif @@ -937,13 +940,13 @@ subroutine read_mesh_databases_hdf5() if (I_should_read_the_database) then call h5_read_dataset_collect_hyperslab("offset_num_abs_boundary_faces",offset_num_abs_boundary_faces, (/0/), H5_COL) endif - call bcast_all_i_array_for_database(offset_num_abs_boundary_faces, size(offset_num_abs_boundary_faces)) + call bcast_all_i_array_for_database(offset_num_abs_boundary_faces, size(offset_num_abs_boundary_faces,kind=4)) if (sum(offset_num_abs_boundary_faces) > 0) then if (I_should_read_the_database) then call h5_read_dataset_collect_hyperslab("offset_nglob_xy",offset_nglob_xy, (/0/), H5_COL) endif - call bcast_all_i_array_for_database(offset_nglob_xy, size(offset_nglob_xy)) + call bcast_all_i_array_for_database(offset_nglob_xy, size(offset_nglob_xy,kind=4)) if (I_should_read_the_database) then dsetname = "abs_boundary_ispec" @@ -960,13 +963,13 @@ subroutine read_mesh_databases_hdf5() (/0,0,sum(offset_num_abs_boundary_faces(0:myrank-1))/), H5_COL) endif if (size(abs_boundary_ispec) > 0) & - call bcast_all_i_for_database(abs_boundary_ispec(1), size(abs_boundary_ispec)) + call bcast_all_i_for_database(abs_boundary_ispec(1), size(abs_boundary_ispec,kind=4)) if (size(abs_boundary_ijk) > 0) & - call bcast_all_i_for_database(abs_boundary_ijk(1,1,1), size(abs_boundary_ijk)) + call bcast_all_i_for_database(abs_boundary_ijk(1,1,1), size(abs_boundary_ijk,kind=4)) if (size(abs_boundary_jacobian2Dw) > 0) & - call bcast_all_cr_for_database(abs_boundary_jacobian2Dw(1,1), size(abs_boundary_jacobian2Dw)) + call bcast_all_cr_for_database(abs_boundary_jacobian2Dw(1,1), size(abs_boundary_jacobian2Dw,kind=4)) if (size(abs_boundary_normal) > 0) & - call bcast_all_cr_for_database(abs_boundary_normal(1,1,1), size(abs_boundary_normal)) + call bcast_all_cr_for_database(abs_boundary_normal(1,1,1), size(abs_boundary_normal,kind=4)) if (STACEY_ABSORBING_CONDITIONS .and. (.not. PML_CONDITIONS)) then ! store mass matrix contributions @@ -982,9 +985,9 @@ subroutine read_mesh_databases_hdf5() call h5_read_dataset_collect_hyperslab(dsetname, rmassz(1:offset_nglob_xy(myrank)), & (/sum(offset_nglob_xy(0:myrank-1))/),H5_COL) endif - if (size(rmassx) > 0) call bcast_all_cr_for_database(rmassx(1), size(rmassx)) - if (size(rmassy) > 0) call bcast_all_cr_for_database(rmassy(1), size(rmassy)) - if (size(rmassz) > 0) call bcast_all_cr_for_database(rmassz(1), size(rmassz)) + if (size(rmassx) > 0) call bcast_all_cr_for_database(rmassx(1), size(rmassx,kind=4)) + if (size(rmassy) > 0) call bcast_all_cr_for_database(rmassy(1), size(rmassy,kind=4)) + if (size(rmassz) > 0) call bcast_all_cr_for_database(rmassz(1), size(rmassz,kind=4)) endif if (ACOUSTIC_SIMULATION) then if (I_should_read_the_database) then @@ -992,7 +995,7 @@ subroutine read_mesh_databases_hdf5() call h5_read_dataset_collect_hyperslab(dsetname, rmassz_acoustic(1:offset_nglob_xy(myrank)), & (/sum(offset_nglob_xy(0:myrank-1))/),H5_COL) endif - if (size(rmassz_acoustic) > 0) call bcast_all_cr_for_database(rmassz_acoustic(1), size(rmassz_acoustic)) + if (size(rmassz_acoustic) > 0) call bcast_all_cr_for_database(rmassz_acoustic(1), size(rmassz_acoustic,kind=4)) endif endif endif @@ -1062,12 +1065,12 @@ subroutine read_mesh_databases_hdf5() call h5_read_dataset_collect_hyperslab(dsetname, ibelm_top, (/sum(offset_nspec2D_top_ext(0:myrank-1))/),H5_COL) endif endif - if (size(ibelm_xmin) > 0) call bcast_all_i_for_database(ibelm_xmin(1), size(ibelm_xmin)) - if (size(ibelm_xmax) > 0) call bcast_all_i_for_database(ibelm_xmax(1), size(ibelm_xmax)) - if (size(ibelm_ymin) > 0) call bcast_all_i_for_database(ibelm_ymin(1), size(ibelm_ymin)) - if (size(ibelm_ymax) > 0) call bcast_all_i_for_database(ibelm_ymax(1), size(ibelm_ymax)) - if (size(ibelm_bottom) > 0) call bcast_all_i_for_database(ibelm_bottom(1), size(ibelm_bottom)) - if (size(ibelm_top) > 0) call bcast_all_i_for_database(ibelm_top(1), size(ibelm_top)) + if (size(ibelm_xmin) > 0) call bcast_all_i_for_database(ibelm_xmin(1), size(ibelm_xmin,kind=4)) + if (size(ibelm_xmax) > 0) call bcast_all_i_for_database(ibelm_xmax(1), size(ibelm_xmax,kind=4)) + if (size(ibelm_ymin) > 0) call bcast_all_i_for_database(ibelm_ymin(1), size(ibelm_ymin,kind=4)) + if (size(ibelm_ymax) > 0) call bcast_all_i_for_database(ibelm_ymax(1), size(ibelm_ymax,kind=4)) + if (size(ibelm_bottom) > 0) call bcast_all_i_for_database(ibelm_bottom(1), size(ibelm_bottom,kind=4)) + if (size(ibelm_top) > 0) call bcast_all_i_for_database(ibelm_top(1), size(ibelm_top,kind=4)) ! free surface if (I_should_read_the_database) then @@ -1091,7 +1094,7 @@ subroutine read_mesh_databases_hdf5() if (I_should_read_the_database) then call h5_read_dataset_collect_hyperslab("offset_num_free_surface_faces",offset_num_free_surface_faces, (/0/), H5_COL) endif - call bcast_all_i_array_for_database(offset_num_free_surface_faces, size(offset_num_free_surface_faces)) + call bcast_all_i_array_for_database(offset_num_free_surface_faces, size(offset_num_free_surface_faces,kind=4)) if (sum(offset_num_free_surface_faces) > 0) then if (I_should_read_the_database) then @@ -1109,13 +1112,13 @@ subroutine read_mesh_databases_hdf5() (/0,0,sum(offset_num_free_surface_faces(0:myrank-1))/),H5_COL) endif if (size(free_surface_ispec) > 0) & - call bcast_all_i_for_database(free_surface_ispec(1), size(free_surface_ispec)) + call bcast_all_i_for_database(free_surface_ispec(1), size(free_surface_ispec,kind=4)) if (size(free_surface_ijk) > 0) & - call bcast_all_i_for_database(free_surface_ijk(1,1,1), size(free_surface_ijk)) + call bcast_all_i_for_database(free_surface_ijk(1,1,1), size(free_surface_ijk,kind=4)) if (size(free_surface_jacobian2Dw) > 0) & - call bcast_all_cr_for_database(free_surface_jacobian2Dw(1,1), size(free_surface_jacobian2Dw)) + call bcast_all_cr_for_database(free_surface_jacobian2Dw(1,1), size(free_surface_jacobian2Dw,kind=4)) if (size(free_surface_normal) > 0) & - call bcast_all_cr_for_database(free_surface_normal(1,1,1), size(free_surface_normal)) + call bcast_all_cr_for_database(free_surface_normal(1,1,1), size(free_surface_normal,kind=4)) endif ! acoustic-elastic coupling surface @@ -1140,7 +1143,7 @@ subroutine read_mesh_databases_hdf5() if (I_should_read_the_database) then call h5_read_dataset_collect_hyperslab("offset_num_coupling_ac_el_faces",offset_num_coupling_ac_el_faces,(/0/),H5_COL) endif - call bcast_all_i_array_for_database(offset_num_coupling_ac_el_faces, size(offset_num_coupling_ac_el_faces)) + call bcast_all_i_array_for_database(offset_num_coupling_ac_el_faces, size(offset_num_coupling_ac_el_faces,kind=4)) if (sum(offset_num_coupling_ac_el_faces) > 0) then if (I_should_read_the_database) then @@ -1158,13 +1161,13 @@ subroutine read_mesh_databases_hdf5() (/0,0,sum(offset_num_coupling_ac_el_faces(0:myrank-1))/),H5_COL) endif if (size(coupling_ac_el_ispec) > 0) & - call bcast_all_i_for_database(coupling_ac_el_ispec(1), size(coupling_ac_el_ispec)) + call bcast_all_i_for_database(coupling_ac_el_ispec(1), size(coupling_ac_el_ispec,kind=4)) if (size(coupling_ac_el_ijk) > 0) & - call bcast_all_i_for_database(coupling_ac_el_ijk(1,1,1), size(coupling_ac_el_ijk)) + call bcast_all_i_for_database(coupling_ac_el_ijk(1,1,1), size(coupling_ac_el_ijk,kind=4)) if (size(coupling_ac_el_jacobian2Dw) > 0) & - call bcast_all_cr_for_database(coupling_ac_el_jacobian2Dw(1,1), size(coupling_ac_el_jacobian2Dw)) + call bcast_all_cr_for_database(coupling_ac_el_jacobian2Dw(1,1), size(coupling_ac_el_jacobian2Dw,kind=4)) if (size(coupling_ac_el_normal) > 0) & - call bcast_all_cr_for_database(coupling_ac_el_normal(1,1,1), size(coupling_ac_el_normal)) + call bcast_all_cr_for_database(coupling_ac_el_normal(1,1,1), size(coupling_ac_el_normal,kind=4)) endif ! acoustic-poroelastic coupling surface @@ -1189,7 +1192,7 @@ subroutine read_mesh_databases_hdf5() if (I_should_read_the_database) then call h5_read_dataset_collect_hyperslab("offset_num_coupling_ac_po_faces",offset_num_coupling_ac_po_faces, (/0/), H5_COL) endif - call bcast_all_i_array_for_database(offset_num_coupling_ac_po_faces,size(offset_num_coupling_ac_po_faces)) + call bcast_all_i_array_for_database(offset_num_coupling_ac_po_faces,size(offset_num_coupling_ac_po_faces,kind=4)) if (sum(offset_num_coupling_ac_po_faces) > 0) then if (I_should_read_the_database) then @@ -1207,13 +1210,13 @@ subroutine read_mesh_databases_hdf5() (/0,0,sum(offset_num_coupling_ac_po_faces(0:myrank-1))/),H5_COL) endif if (size(coupling_ac_po_ispec) > 0) & - call bcast_all_i_for_database(coupling_ac_po_ispec(1), size(coupling_ac_po_ispec)) + call bcast_all_i_for_database(coupling_ac_po_ispec(1), size(coupling_ac_po_ispec,kind=4)) if (size(coupling_ac_po_ijk) > 0) & - call bcast_all_i_for_database(coupling_ac_po_ijk(1,1,1), size(coupling_ac_po_ijk)) + call bcast_all_i_for_database(coupling_ac_po_ijk(1,1,1), size(coupling_ac_po_ijk,kind=4)) if (size(coupling_ac_po_jacobian2Dw) > 0) & - call bcast_all_cr_for_database(coupling_ac_po_jacobian2Dw(1,1), size(coupling_ac_po_jacobian2Dw)) + call bcast_all_cr_for_database(coupling_ac_po_jacobian2Dw(1,1), size(coupling_ac_po_jacobian2Dw,kind=4)) if (size(coupling_ac_po_normal) > 0) & - call bcast_all_cr_for_database(coupling_ac_po_normal(1,1,1), size(coupling_ac_po_normal)) + call bcast_all_cr_for_database(coupling_ac_po_normal(1,1,1), size(coupling_ac_po_normal,kind=4)) endif ! elastic-poroelastic coupling surface @@ -1243,7 +1246,7 @@ subroutine read_mesh_databases_hdf5() if (I_should_read_the_database) then call h5_read_dataset_collect_hyperslab("offset_num_coupling_el_po_faces",offset_num_coupling_el_po_faces, (/0/), H5_COL) endif - call bcast_all_i_array_for_database(offset_num_coupling_el_po_faces,size(offset_num_coupling_el_po_faces)) + call bcast_all_i_array_for_database(offset_num_coupling_el_po_faces,size(offset_num_coupling_el_po_faces,kind=4)) if (sum(offset_num_coupling_el_po_faces) > 0) then if (I_should_read_the_database) then @@ -1267,17 +1270,17 @@ subroutine read_mesh_databases_hdf5() (/0,0,sum(offset_num_coupling_el_po_faces(0:myrank-1))/),H5_COL) endif if (size(coupling_el_po_ispec) > 0) & - call bcast_all_i_for_database(coupling_el_po_ispec(1), size(coupling_el_po_ispec)) + call bcast_all_i_for_database(coupling_el_po_ispec(1), size(coupling_el_po_ispec,kind=4)) if (size(coupling_po_el_ispec) > 0) & - call bcast_all_i_for_database(coupling_po_el_ispec(1), size(coupling_po_el_ispec)) + call bcast_all_i_for_database(coupling_po_el_ispec(1), size(coupling_po_el_ispec,kind=4)) if (size(coupling_el_po_ijk) > 0) & - call bcast_all_i_for_database(coupling_el_po_ijk(1,1,1), size(coupling_el_po_ijk)) + call bcast_all_i_for_database(coupling_el_po_ijk(1,1,1), size(coupling_el_po_ijk,kind=4)) if (size(coupling_po_el_ijk) > 0) & - call bcast_all_i_for_database(coupling_po_el_ijk(1,1,1), size(coupling_po_el_ijk)) + call bcast_all_i_for_database(coupling_po_el_ijk(1,1,1), size(coupling_po_el_ijk,kind=4)) if (size(coupling_el_po_jacobian2Dw) > 0) & - call bcast_all_cr_for_database(coupling_el_po_jacobian2Dw(1,1), size(coupling_el_po_jacobian2Dw)) + call bcast_all_cr_for_database(coupling_el_po_jacobian2Dw(1,1), size(coupling_el_po_jacobian2Dw,kind=4)) if (size(coupling_el_po_normal) > 0) & - call bcast_all_cr_for_database(coupling_el_po_normal(1,1,1), size(coupling_el_po_normal)) + call bcast_all_cr_for_database(coupling_el_po_normal(1,1,1), size(coupling_el_po_normal,kind=4)) endif ! MPI interfaces @@ -1297,7 +1300,7 @@ subroutine read_mesh_databases_hdf5() if (I_should_read_the_database) then call h5_read_dataset_collect_hyperslab("offset_num_interfaces_ext_mesh",offset_num_interfaces_ext_mesh, (/0/), H5_COL) endif - call bcast_all_i_array_for_database(offset_num_interfaces_ext_mesh,size(offset_num_interfaces_ext_mesh)) + call bcast_all_i_array_for_database(offset_num_interfaces_ext_mesh,size(offset_num_interfaces_ext_mesh,kind=4)) if (sum(offset_num_interfaces_ext_mesh) > 0) then if (I_should_read_the_database) then @@ -1322,11 +1325,11 @@ subroutine read_mesh_databases_hdf5() (/0,sum(offset_num_interfaces_ext_mesh(0:myrank-1))/),H5_COL) endif if (size(my_neighbors_ext_mesh) > 0) & - call bcast_all_i_for_database(my_neighbors_ext_mesh(1), size(my_neighbors_ext_mesh)) + call bcast_all_i_for_database(my_neighbors_ext_mesh(1), size(my_neighbors_ext_mesh,kind=4)) if (size(nibool_interfaces_ext_mesh) > 0) & - call bcast_all_i_for_database(nibool_interfaces_ext_mesh(1), size(nibool_interfaces_ext_mesh)) + call bcast_all_i_for_database(nibool_interfaces_ext_mesh(1), size(nibool_interfaces_ext_mesh,kind=4)) if (size(ibool_interfaces_ext_mesh) > 0) & - call bcast_all_i_for_database(ibool_interfaces_ext_mesh(1,1), size(ibool_interfaces_ext_mesh)) + call bcast_all_i_for_database(ibool_interfaces_ext_mesh(1,1), size(ibool_interfaces_ext_mesh,kind=4)) else ! no interfaces max_nibool_interfaces_ext_mesh = 0 @@ -1382,27 +1385,27 @@ subroutine read_mesh_databases_hdf5() dsetname = "c66store" call h5_read_dataset_collect_hyperslab(dsetname, c66store, (/0,0,0,sum(offset_nspec_aniso(0:myrank-1))/), H5_COL) endif - if (size(c11store) > 0) call bcast_all_cr_for_database(c11store(1,1,1,1), size(c11store)) - if (size(c12store) > 0) call bcast_all_cr_for_database(c12store(1,1,1,1), size(c12store)) - if (size(c13store) > 0) call bcast_all_cr_for_database(c13store(1,1,1,1), size(c13store)) - if (size(c14store) > 0) call bcast_all_cr_for_database(c14store(1,1,1,1), size(c14store)) - if (size(c15store) > 0) call bcast_all_cr_for_database(c15store(1,1,1,1), size(c15store)) - if (size(c16store) > 0) call bcast_all_cr_for_database(c16store(1,1,1,1), size(c16store)) - if (size(c22store) > 0) call bcast_all_cr_for_database(c22store(1,1,1,1), size(c22store)) - if (size(c23store) > 0) call bcast_all_cr_for_database(c23store(1,1,1,1), size(c23store)) - if (size(c24store) > 0) call bcast_all_cr_for_database(c24store(1,1,1,1), size(c24store)) - if (size(c25store) > 0) call bcast_all_cr_for_database(c25store(1,1,1,1), size(c25store)) - if (size(c26store) > 0) call bcast_all_cr_for_database(c26store(1,1,1,1), size(c26store)) - if (size(c33store) > 0) call bcast_all_cr_for_database(c33store(1,1,1,1), size(c33store)) - if (size(c34store) > 0) call bcast_all_cr_for_database(c34store(1,1,1,1), size(c34store)) - if (size(c35store) > 0) call bcast_all_cr_for_database(c35store(1,1,1,1), size(c35store)) - if (size(c36store) > 0) call bcast_all_cr_for_database(c36store(1,1,1,1), size(c36store)) - if (size(c44store) > 0) call bcast_all_cr_for_database(c44store(1,1,1,1), size(c44store)) - if (size(c45store) > 0) call bcast_all_cr_for_database(c45store(1,1,1,1), size(c45store)) - if (size(c46store) > 0) call bcast_all_cr_for_database(c46store(1,1,1,1), size(c46store)) - if (size(c55store) > 0) call bcast_all_cr_for_database(c55store(1,1,1,1), size(c55store)) - if (size(c56store) > 0) call bcast_all_cr_for_database(c56store(1,1,1,1), size(c56store)) - if (size(c66store) > 0) call bcast_all_cr_for_database(c66store(1,1,1,1), size(c66store)) + if (size(c11store) > 0) call bcast_all_cr_for_database(c11store(1,1,1,1), size(c11store,kind=4)) + if (size(c12store) > 0) call bcast_all_cr_for_database(c12store(1,1,1,1), size(c12store,kind=4)) + if (size(c13store) > 0) call bcast_all_cr_for_database(c13store(1,1,1,1), size(c13store,kind=4)) + if (size(c14store) > 0) call bcast_all_cr_for_database(c14store(1,1,1,1), size(c14store,kind=4)) + if (size(c15store) > 0) call bcast_all_cr_for_database(c15store(1,1,1,1), size(c15store,kind=4)) + if (size(c16store) > 0) call bcast_all_cr_for_database(c16store(1,1,1,1), size(c16store,kind=4)) + if (size(c22store) > 0) call bcast_all_cr_for_database(c22store(1,1,1,1), size(c22store,kind=4)) + if (size(c23store) > 0) call bcast_all_cr_for_database(c23store(1,1,1,1), size(c23store,kind=4)) + if (size(c24store) > 0) call bcast_all_cr_for_database(c24store(1,1,1,1), size(c24store,kind=4)) + if (size(c25store) > 0) call bcast_all_cr_for_database(c25store(1,1,1,1), size(c25store,kind=4)) + if (size(c26store) > 0) call bcast_all_cr_for_database(c26store(1,1,1,1), size(c26store,kind=4)) + if (size(c33store) > 0) call bcast_all_cr_for_database(c33store(1,1,1,1), size(c33store,kind=4)) + if (size(c34store) > 0) call bcast_all_cr_for_database(c34store(1,1,1,1), size(c34store,kind=4)) + if (size(c35store) > 0) call bcast_all_cr_for_database(c35store(1,1,1,1), size(c35store,kind=4)) + if (size(c36store) > 0) call bcast_all_cr_for_database(c36store(1,1,1,1), size(c36store,kind=4)) + if (size(c44store) > 0) call bcast_all_cr_for_database(c44store(1,1,1,1), size(c44store,kind=4)) + if (size(c45store) > 0) call bcast_all_cr_for_database(c45store(1,1,1,1), size(c45store,kind=4)) + if (size(c46store) > 0) call bcast_all_cr_for_database(c46store(1,1,1,1), size(c46store,kind=4)) + if (size(c55store) > 0) call bcast_all_cr_for_database(c55store(1,1,1,1), size(c55store,kind=4)) + if (size(c56store) > 0) call bcast_all_cr_for_database(c56store(1,1,1,1), size(c56store,kind=4)) + if (size(c66store) > 0) call bcast_all_cr_for_database(c66store(1,1,1,1), size(c66store,kind=4)) endif ! inner / outer elements @@ -1415,7 +1418,7 @@ subroutine read_mesh_databases_hdf5() dsetname = "ispec_is_inner" call h5_read_dataset_collect_hyperslab(dsetname, ispec_is_inner,(/sum(offset_nspec(0:myrank-1))/),H5_COL) endif - if (size(ispec_is_inner) > 0) call bcast_all_l_for_database(ispec_is_inner(1), size(ispec_is_inner)) + if (size(ispec_is_inner) > 0) call bcast_all_l_for_database(ispec_is_inner(1), size(ispec_is_inner,kind=4)) if (ACOUSTIC_SIMULATION) then if (I_should_read_the_database) then @@ -1445,7 +1448,7 @@ subroutine read_mesh_databases_hdf5() (/sum(offset_num_phase_ispec_acoustic(0:myrank-1)),0/),H5_COL) endif if (size(phase_ispec_inner_acoustic) > 0) & - call bcast_all_i_for_database(phase_ispec_inner_acoustic(1,1), size(phase_ispec_inner_acoustic)) + call bcast_all_i_for_database(phase_ispec_inner_acoustic(1,1), size(phase_ispec_inner_acoustic,kind=4)) endif endif @@ -1471,7 +1474,7 @@ subroutine read_mesh_databases_hdf5() if (I_should_read_the_database) then call h5_read_dataset_collect_hyperslab("offset_num_phase_ispec_elastic",offset_num_phase_ispec_elastic, (/0/), H5_COL) endif - call bcast_all_i_array_for_database(offset_num_phase_ispec_elastic,size(offset_num_phase_ispec_elastic)) + call bcast_all_i_array_for_database(offset_num_phase_ispec_elastic,size(offset_num_phase_ispec_elastic,kind=4)) if (sum(offset_num_phase_ispec_elastic) > 0) then if (I_should_read_the_database) then @@ -1480,7 +1483,7 @@ subroutine read_mesh_databases_hdf5() (/sum(offset_num_phase_ispec_elastic(0:myrank-1)),0/),H5_COL) endif if (size(phase_ispec_inner_elastic) > 0) & - call bcast_all_i_for_database(phase_ispec_inner_elastic(1,1), size(phase_ispec_inner_elastic)) + call bcast_all_i_for_database(phase_ispec_inner_elastic(1,1), size(phase_ispec_inner_elastic,kind=4)) endif endif @@ -1507,7 +1510,7 @@ subroutine read_mesh_databases_hdf5() call h5_read_dataset_collect_hyperslab("offset_num_phase_ispec_poroelastic", & offset_num_phase_ispec_poroelastic, (/0/), H5_COL) endif - call bcast_all_i_array_for_database(offset_num_phase_ispec_poroelastic,size(offset_num_phase_ispec_poroelastic)) + call bcast_all_i_array_for_database(offset_num_phase_ispec_poroelastic,size(offset_num_phase_ispec_poroelastic,kind=4)) if (sum(offset_num_phase_ispec_poroelastic) > 0) then if (I_should_read_the_database) then dsetname = "phase_ispec_inner_poroelastic" @@ -1515,7 +1518,7 @@ subroutine read_mesh_databases_hdf5() (/sum(offset_num_phase_ispec_poroelastic(0:myrank-1)),0/),H5_COL) endif if (size(phase_ispec_inner_poroelastic) > 0) & - call bcast_all_i_for_database(phase_ispec_inner_poroelastic(1,1), size(phase_ispec_inner_poroelastic)) + call bcast_all_i_for_database(phase_ispec_inner_poroelastic(1,1), size(phase_ispec_inner_poroelastic,kind=4)) endif endif @@ -1545,7 +1548,7 @@ subroutine read_mesh_databases_hdf5() num_elem_colors_acoustic,(/sum(offset_num_colors_acoustic(0:myrank-1))/),H5_COL) endif if (size(num_elem_colors_acoustic) > 0) & - call bcast_all_i_for_database(num_elem_colors_acoustic(1), size(num_elem_colors_acoustic)) + call bcast_all_i_for_database(num_elem_colors_acoustic(1), size(num_elem_colors_acoustic,kind=4)) endif ! elastic domain colors if (ELASTIC_SIMULATION) then @@ -1571,7 +1574,7 @@ subroutine read_mesh_databases_hdf5() (/sum(offset_num_colors_elastic(0:myrank-1))/),H5_COL) endif if (size(num_elem_colors_elastic) > 0) & - call bcast_all_i_for_database(num_elem_colors_elastic(1), size(num_elem_colors_elastic)) + call bcast_all_i_for_database(num_elem_colors_elastic(1), size(num_elem_colors_elastic,kind=4)) endif else ! allocates dummy arrays @@ -1614,8 +1617,8 @@ subroutine read_mesh_databases_hdf5() (/sum(offset_nglob_ab(0:myrank-1))/), H5_COL) endif call bcast_all_i_for_database(nfaces_surface, 1) - call bcast_all_l_for_database(ispec_is_surface_external_mesh(1), size(ispec_is_surface_external_mesh)) - call bcast_all_l_for_database(iglob_is_surface_external_mesh(1), size(iglob_is_surface_external_mesh)) + call bcast_all_l_for_database(ispec_is_surface_external_mesh(1), size(ispec_is_surface_external_mesh,kind=4)) + call bcast_all_l_for_database(iglob_is_surface_external_mesh(1), size(iglob_is_surface_external_mesh,kind=4)) ! for mesh adjacency if (I_should_read_the_database) then @@ -1640,8 +1643,8 @@ subroutine read_mesh_databases_hdf5() call h5_read_dataset_collect_hyperslab(dsetname, neighbors_adjncy, & (/sum_offset_neighbors_adjncy_this_proc/), H5_COL) endif - call bcast_all_i_for_database(neighbors_xadj(1), size(neighbors_xadj)) - call bcast_all_i_for_database(neighbors_adjncy(1), size(neighbors_adjncy)) + call bcast_all_i_for_database(neighbors_xadj(1), size(neighbors_xadj,kind=4)) + call bcast_all_i_for_database(neighbors_adjncy(1), size(neighbors_adjncy,kind=4)) ! done reading database if (I_should_read_the_database) then diff --git a/src/specfem3D/save_forward_arrays.f90 b/src/specfem3D/save_forward_arrays.f90 index a4f550409..f5b1811c5 100644 --- a/src/specfem3D/save_forward_arrays.f90 +++ b/src/specfem3D/save_forward_arrays.f90 @@ -127,7 +127,7 @@ subroutine save_forward_arrays_undoatt() if (ATTENUATION) then ! only memory variables needed call transfer_rmemory_from_device(Mesh_pointer,R_xx,R_yy,R_xy,R_xz,R_yz, & - R_trace,size(R_xx)) + R_trace,size(R_xx,kind=4)) endif endif endif diff --git a/src/specfem3D/vtk_window.F90 b/src/specfem3D/vtk_window.F90 index d2ef2f7a4..fd68dda84 100644 --- a/src/specfem3D/vtk_window.F90 +++ b/src/specfem3D/vtk_window.F90 @@ -1040,12 +1040,12 @@ subroutine vtk_window_update() if (NPROC > 1) then if (myrank == 0) then ! gather data - call gatherv_all_cr(vtkdata,size(vtkdata), & + call gatherv_all_cr(vtkdata,size(vtkdata,kind=4), & vtkdata_all,vtkdata_points_all,vtkdata_offset_all, & vtkdata_numpoints_all,NPROC) else ! all other process just send data - call gatherv_all_cr(vtkdata,size(vtkdata), & + call gatherv_all_cr(vtkdata,size(vtkdata,kind=4), & dummy,vtkdata_points_all,vtkdata_offset_all, & 1,NPROC) endif diff --git a/src/specfem3D/write_movie_output_HDF5.F90 b/src/specfem3D/write_movie_output_HDF5.F90 index 04d359f02..31e176737 100644 --- a/src/specfem3D/write_movie_output_HDF5.F90 +++ b/src/specfem3D/write_movie_output_HDF5.F90 @@ -201,7 +201,7 @@ subroutine wmo_save_shakemap_hdf5() fname_h5_data_shake = trim(OUTPUT_FILES) // "/shakemap.h5" ! get the offset info from main rank - call bcast_all_i(faces_surface_offset,size(faces_surface_offset)) + call bcast_all_i(faces_surface_offset,size(faces_surface_offset,kind=4)) if (myrank == 0) then size_surf_array = size(store_val_x_all) @@ -474,7 +474,7 @@ subroutine wmo_movie_save_surface_snapshot_hdf5() ! get the offset info from main rank if (it == NTSTEP_BETWEEN_FRAMES) then - call bcast_all_i(faces_surface_offset,size(faces_surface_offset)) + call bcast_all_i(faces_surface_offset,size(faces_surface_offset,kind=4)) endif if (myrank == 0) then From f1aed570a1cb9935c2461bbe2205d4470a9f5abc Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Wed, 30 Oct 2024 16:00:37 +0100 Subject: [PATCH 5/5] updates receiver check --- .../setup_mesh_adjacency.f90 | 3 +- src/specfem3D/locate_receivers.f90 | 53 ++++++++++++------- 2 files changed, 37 insertions(+), 19 deletions(-) diff --git a/src/generate_databases/setup_mesh_adjacency.f90 b/src/generate_databases/setup_mesh_adjacency.f90 index 23d2a1ecd..63c5af864 100644 --- a/src/generate_databases/setup_mesh_adjacency.f90 +++ b/src/generate_databases/setup_mesh_adjacency.f90 @@ -255,7 +255,8 @@ subroutine setup_mesh_adjacency() else ! no neighbors ! warning - print *,'*** Warning: found mesh element with no neighbors : slice ',myrank,' - element ',ispec_ref,' ***' + print *,'*** Warning: found mesh element with no neighbors : slice ',myrank, & + ' - element ',ispec_ref,'out of',NSPEC_AB,' ***' endif ! again loop to get neighbors of neighbors diff --git a/src/specfem3D/locate_receivers.f90 b/src/specfem3D/locate_receivers.f90 index b3301cbdf..fffc050bc 100644 --- a/src/specfem3D/locate_receivers.f90 +++ b/src/specfem3D/locate_receivers.f90 @@ -39,7 +39,8 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) xi_receiver,eta_receiver,gamma_receiver,nu_rec, & station_name,network_name, & stlat,stlon,stbur, & - x_target_station,y_target_station,z_target_station + x_target_station,y_target_station,z_target_station, & + NPROC ! PML use pml_par, only: is_CPML @@ -55,7 +56,7 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) double precision, allocatable, dimension(:) :: x_target,y_target,z_target double precision, allocatable, dimension(:) :: x_found,y_found,z_found - integer :: irec,ier,ispec + integer :: irec,ier,ispec,islice ! timer MPI double precision, external :: wtime @@ -283,35 +284,51 @@ subroutine locate_receivers(rec_filename,utm_x_source,utm_y_source) call bcast_all_dp(nu_rec,NDIM*NDIM*nrec) call bcast_all_dp(final_distance,nrec) - ! checks if valid receiver element - if (myrank == 0) then - ! locate point might return a zero ispec if point outside/above mesh - do irec = 1,nrec + ! checks if we got valid receiver elements + ! locate point might return a zero ispec if point outside/above mesh + do irec = 1,nrec + ! slice the receiver is in + islice = islice_selected_rec(irec) + + ! checks slice + if (islice < 0 .or. islice > NPROC-1) then + print *,'Error locating station # ',irec,' ',trim(network_name(irec)),' ',trim(station_name(irec)) + print *,' found in an invalid slice: ',islice + call exit_MPI(myrank,'Error something is wrong with the slice number of receiver') + endif + + ! checks found element + if (myrank == islice_selected_rec(irec)) then + ! element index ispec = ispec_selected_rec(irec) ! checks if valid if (ispec < 1 .or. ispec > NSPEC_AB) then ! invalid element print *,'Error locating station # ',irec,' ',trim(network_name(irec)),' ',trim(station_name(irec)) if (SUPPRESS_UTM_PROJECTION) then - print *,' original x: ',sngl(stutm_x(irec)) - print *,' original y: ',sngl(stutm_y(irec)) + print *,' original x: ',sngl(stutm_x(irec)) + print *,' original y: ',sngl(stutm_y(irec)) else - print *,' original UTM x: ',sngl(stutm_x(irec)) - print *,' original UTM y: ',sngl(stutm_y(irec)) + print *,' original UTM x: ',sngl(stutm_x(irec)) + print *,' original UTM y: ',sngl(stutm_y(irec)) endif if (USE_SOURCES_RECEIVERS_Z) then - print *,' original z: ',sngl(stbur(irec)) + print *,' original z: ',sngl(stbur(irec)) else - print *,' original depth: ',sngl(stbur(irec)),' m' + print *,' original depth: ',sngl(stbur(irec)),' m' endif - print *,' only found invalid element: slice ',islice_selected_rec(irec) - print *,' ispec ',ispec_selected_rec(irec) - print *,' domain ',idomain(irec) - print *,' final_distance: ',final_distance(irec) + print *,' found in an invalid element: slice :',islice_selected_rec(irec) + print *,' ispec :',ispec_selected_rec(irec),'out of ',NSPEC_AB + print *,' domain :',idomain(irec) + print *,' xi/eta/gamma :',xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec) + print *,' final_distance: ',final_distance(irec) + print * + print *,'Please check your receiver position in file DATA/STATIONS, and move it closer the mesh geometry - exiting...' + print * call exit_MPI(myrank,'Error locating receiver') endif - enddo - endif + endif + enddo call synchronize_all() ! warning if receiver in C-PML region