[cig-commits] r13206 - in seismo/3D/SPECFEM3D/trunk: . UTILS
nlegoff at geodynamics.org
nlegoff at geodynamics.org
Fri Oct 31 12:04:49 PDT 2008
Author: nlegoff
Date: 2008-10-31 12:04:49 -0700 (Fri, 31 Oct 2008)
New Revision: 13206
Added:
seismo/3D/SPECFEM3D/trunk/UTILS/external_mesh/
seismo/3D/SPECFEM3D/trunk/compute_forces.f90
seismo/3D/SPECFEM3D/trunk/sort_array_coordinates.f90
Modified:
seismo/3D/SPECFEM3D/trunk/Makefile.in
seismo/3D/SPECFEM3D/trunk/assemble_MPI_scalar.f90
seismo/3D/SPECFEM3D/trunk/assemble_MPI_vector.f90
seismo/3D/SPECFEM3D/trunk/constants.h.in
seismo/3D/SPECFEM3D/trunk/create_movie_AVS_DX.f90
seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/locate_receivers.f90
seismo/3D/SPECFEM3D/trunk/locate_source.f90
seismo/3D/SPECFEM3D/trunk/meshfem3D.f90
seismo/3D/SPECFEM3D/trunk/parallel.f90
seismo/3D/SPECFEM3D/trunk/read_arrays_solver.f90
seismo/3D/SPECFEM3D/trunk/specfem3D.f90
Log:
added use of external meshes. See UTILS/external_mesh for more info.
Modified: seismo/3D/SPECFEM3D/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/Makefile.in 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/Makefile.in 2008-10-31 19:04:49 UTC (rev 13206)
@@ -122,6 +122,7 @@
$O/save_arrays_solver.o \
$O/save_header_file.o \
$O/socal_model.o \
+ $O/sort_array_coordinates.o \
$O/utm_geo.o \
$O/write_AVS_DX_global_data.o \
$O/write_AVS_DX_global_faces_data.o \
@@ -129,6 +130,7 @@
$O/write_AVS_DX_surface_data.o \
$O/write_seismograms.o \
$O/compute_boundary_kernel.o \
+ $O/compute_forces.o \
$(EMPTY_MACRO)
# solver objects with statically allocated arrays; dependent upon
@@ -420,6 +422,9 @@
$O/socal_model.o: constants.h socal_model.f90
${FCCOMPILE_CHECK} -c -o $O/socal_model.o socal_model.f90
+$O/sort_array_coordinates.o: constants.h sort_array_coordinates.f90
+ ${FCCOMPILE_CHECK} -c -o $O/sort_array_coordinates.o sort_array_coordinates.f90
+
$O/aniso_model.o: constants.h aniso_model.f90
${FCCOMPILE_CHECK} -c -o $O/aniso_model.o aniso_model.f90
@@ -480,6 +485,9 @@
$O/compute_boundary_kernel.o: constants.h compute_boundary_kernel.f90
${FCCOMPILE_CHECK} -c -o $O/compute_boundary_kernel.o compute_boundary_kernel.f90
+$O/compute_forces.o: constants.h compute_forces.f90
+ ${FCCOMPILE_CHECK} -c -o $O/compute_forces.o compute_forces.f90
+
$O/combine_vol_data.o: constants.h combine_vol_data.f90
${FCCOMPILE_CHECK} -c -o $O/combine_vol_data.o combine_vol_data.f90
Copied: seismo/3D/SPECFEM3D/trunk/UTILS/external_mesh (from rev 13200, seismo/3D/SPECFEM3D/branches/update_temporary/UTILS/external_mesh)
Modified: seismo/3D/SPECFEM3D/trunk/assemble_MPI_scalar.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/assemble_MPI_scalar.f90 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/assemble_MPI_scalar.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -196,3 +196,172 @@
end subroutine assemble_MPI_scalar
+!
+!----
+!
+
+ subroutine assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,array_val, &
+ buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh &
+ )
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+! array to assemble
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: array_val
+
+ integer :: NPROC
+ integer :: NGLOB_AB
+
+ real(kind=CUSTOM_REAL), dimension(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: &
+ buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh
+
+ integer :: ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbours_ext_mesh
+ integer, dimension(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh
+
+ integer ipoin,iinterface
+ integer sender,receiver
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! here we have to assemble all the contributions between partitions using MPI
+
+! assemble only if more than one partition
+ if(NPROC > 1) then
+
+! partition border copy into the buffer
+ do iinterface = 1, ninterfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ buffer_send_scalar_ext_mesh(ipoin,iinterface) = array_val(ibool_interfaces_ext_mesh(ipoin,iinterface))
+ enddo
+ enddo
+
+! send messages
+ do iinterface = 1, ninterfaces_ext_mesh
+ call issend_cr(buffer_send_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
+ nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_send_scalar_ext_mesh(iinterface) &
+ )
+ call irecv_cr(buffer_recv_scalar_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
+ nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_recv_scalar_ext_mesh(iinterface) &
+ )
+ enddo
+
+! wait for communications completion
+ do iinterface = 1, ninterfaces_ext_mesh
+ call wait_req(request_recv_scalar_ext_mesh(iinterface))
+ enddo
+
+! adding contributions of neighbours
+ do iinterface = 1, ninterfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
+ array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_scalar_ext_mesh(ipoin,iinterface)
+ enddo
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_scalar_ext_mesh
+
+!
+!----
+!
+
+ subroutine assemble_MPI_scalar_i_ext_mesh(NPROC,NGLOB_AB,array_val, &
+ buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh &
+ )
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+! array to assemble
+ integer, dimension(NGLOB_AB) :: array_val
+
+ integer :: NPROC
+ integer :: NGLOB_AB
+
+ integer, dimension(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: &
+ buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh
+
+ integer :: ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbours_ext_mesh
+ integer, dimension(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh
+
+ integer ipoin,iinterface
+ integer sender,receiver
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! here we have to assemble all the contributions between partitions using MPI
+
+! assemble only if more than one partition
+ if(NPROC > 1) then
+
+! partition border copy into the buffer
+ do iinterface = 1, ninterfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ buffer_send_scalar_ext_mesh(ipoin,iinterface) = array_val(ibool_interfaces_ext_mesh(ipoin,iinterface))
+ enddo
+ enddo
+
+! send messages
+ do iinterface = 1, ninterfaces_ext_mesh
+ call issend_i(buffer_send_scalar_ext_mesh(1,iinterface), &
+ nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_send_scalar_ext_mesh(iinterface) &
+ )
+ call irecv_i(buffer_recv_scalar_ext_mesh(1,iinterface), &
+ nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_recv_scalar_ext_mesh(iinterface) &
+ )
+ enddo
+
+! wait for communications completion
+ do iinterface = 1, ninterfaces_ext_mesh
+ call wait_req(request_recv_scalar_ext_mesh(iinterface))
+ enddo
+
+! adding contributions of neighbours
+ do iinterface = 1, ninterfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
+ array_val(ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_scalar_ext_mesh(ipoin,iinterface)
+ enddo
+ enddo
+
+! wait for communications completion (send)
+ do iinterface = 1, ninterfaces_ext_mesh
+ call wait_req(request_send_scalar_ext_mesh(iinterface))
+ enddo
+
+
+ endif
+
+ end subroutine assemble_MPI_scalar_i_ext_mesh
Modified: seismo/3D/SPECFEM3D/trunk/assemble_MPI_vector.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/assemble_MPI_vector.f90 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/assemble_MPI_vector.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -196,3 +196,215 @@
end subroutine assemble_MPI_vector
+!
+!----
+!
+
+ subroutine assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,array_val, &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh &
+ )
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+! array to assemble
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: array_val
+
+ integer :: NPROC
+ integer :: NGLOB_AB
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh
+
+ integer :: ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbours_ext_mesh
+ integer, dimension(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: request_send_vector_ext_mesh,request_recv_vector_ext_mesh
+
+ integer ipoin,iinterface
+ integer sender,receiver
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! here we have to assemble all the contributions between partitions using MPI
+
+! assemble only if more than one partition
+ if(NPROC > 1) then
+
+! partition border copy into the buffer
+ do iinterface = 1, ninterfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ buffer_send_vector_ext_mesh(:,ipoin,iinterface) = array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface))
+ enddo
+ enddo
+
+! send messages
+ do iinterface = 1, ninterfaces_ext_mesh
+ call issend_cr(buffer_send_vector_ext_mesh(1,1,iinterface), &
+ NDIM*nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_send_vector_ext_mesh(iinterface) &
+ )
+ call irecv_cr(buffer_recv_vector_ext_mesh(1,1,iinterface), &
+ NDIM*nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_recv_vector_ext_mesh(iinterface) &
+ )
+ enddo
+
+! wait for communications completion (recv)
+ do iinterface = 1, ninterfaces_ext_mesh
+ call wait_req(request_recv_vector_ext_mesh(iinterface))
+ enddo
+
+! adding contributions of neighbours
+ do iinterface = 1, ninterfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
+ array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_vector_ext_mesh(:,ipoin,iinterface)
+ enddo
+ enddo
+
+! wait for communications completion (send)
+ do iinterface = 1, ninterfaces_ext_mesh
+ call wait_req(request_send_vector_ext_mesh(iinterface))
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_vector_ext_mesh
+
+ subroutine assemble_MPI_vector_ext_mesh_s(NPROC,NGLOB_AB,array_val, &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh &
+ )
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+! array to assemble
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: array_val
+
+ integer :: NPROC
+ integer :: NGLOB_AB
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh
+
+ integer :: ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbours_ext_mesh
+ integer, dimension(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: request_send_vector_ext_mesh,request_recv_vector_ext_mesh
+
+ integer ipoin,iinterface
+ integer sender,receiver
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! here we have to assemble all the contributions between partitions using MPI
+
+! assemble only if more than one partition
+ if(NPROC > 1) then
+
+! partition border copy into the buffer
+ do iinterface = 1, ninterfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ buffer_send_vector_ext_mesh(:,ipoin,iinterface) = array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface))
+ enddo
+ enddo
+
+! send messages
+ do iinterface = 1, ninterfaces_ext_mesh
+ call issend_cr(buffer_send_vector_ext_mesh(1,1,iinterface), &
+ NDIM*nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_send_vector_ext_mesh(iinterface) &
+ )
+ call irecv_cr(buffer_recv_vector_ext_mesh(1,1,iinterface), &
+ NDIM*nibool_interfaces_ext_mesh(iinterface), &
+ my_neighbours_ext_mesh(iinterface), &
+ itag, &
+ request_recv_vector_ext_mesh(iinterface) &
+ )
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_vector_ext_mesh_s
+
+
+ subroutine assemble_MPI_vector_ext_mesh_w(NPROC,NGLOB_AB,array_val, &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh &
+ )
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+! array to assemble
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: array_val
+
+ integer :: NPROC
+ integer :: NGLOB_AB
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh
+
+ integer :: ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbours_ext_mesh
+ integer, dimension(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(ninterfaces_ext_mesh) :: request_send_vector_ext_mesh,request_recv_vector_ext_mesh
+
+ integer ipoin,iinterface
+ integer sender,receiver
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! here we have to assemble all the contributions between partitions using MPI
+
+! assemble only if more than one partition
+ if(NPROC > 1) then
+
+! wait for communications completion (recv)
+ do iinterface = 1, ninterfaces_ext_mesh
+ call wait_req(request_recv_vector_ext_mesh(iinterface))
+ enddo
+
+! adding contributions of neighbours
+ do iinterface = 1, ninterfaces_ext_mesh
+ do ipoin = 1, nibool_interfaces_ext_mesh(iinterface)
+ array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) = &
+ array_val(:,ibool_interfaces_ext_mesh(ipoin,iinterface)) + buffer_recv_vector_ext_mesh(:,ipoin,iinterface)
+ enddo
+ enddo
+
+! wait for communications completion (send)
+ do iinterface = 1, ninterfaces_ext_mesh
+ call wait_req(request_send_vector_ext_mesh(iinterface))
+ enddo
+
+ endif
+
+ end subroutine assemble_MPI_vector_ext_mesh_w
Copied: seismo/3D/SPECFEM3D/trunk/compute_forces.f90 (from rev 13200, seismo/3D/SPECFEM3D/branches/update_temporary/compute_forces.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_forces.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/compute_forces.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -0,0 +1,301 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+subroutine compute_forces(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ kappastore,mustore,jacobian,ibool,ispec_is_inner,phase_is_inner, &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt &
+ )
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ kappastore,mustore,jacobian
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! source
+ integer :: NSOURCES,myrank,it
+ integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
+ double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
+ double precision, dimension(3,3,NSOURCES) :: nu_source
+ double precision, dimension(NSOURCES) :: hdur
+ double precision :: dt
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+
+ integer ispec,iglob
+ integer i,j,k,l
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+ real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+
+ real(kind=CUSTOM_REAL) hp1,hp2,hp3
+ real(kind=CUSTOM_REAL) fac1,fac2,fac3
+
+ real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
+ real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
+ real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+
+ real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
+ real(kind=CUSTOM_REAL) kappal
+
+ integer :: isource
+ double precision :: t0,f0
+
+
+ do ispec = 1,NSPEC_AB
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0.
+ tempx2l = 0.
+ tempx3l = 0.
+
+ tempy1l = 0.
+ tempy2l = 0.
+ tempy3l = 0.
+
+ tempz1l = 0.
+ tempz2l = 0.
+ tempz3l = 0.
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec)
+ tempx1l = tempx1l + displ(1,iglob)*hp1
+ tempy1l = tempy1l + displ(2,iglob)*hp1
+ tempz1l = tempz1l + displ(3,iglob)*hp1
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec)
+ tempx2l = tempx2l + displ(1,iglob)*hp2
+ tempy2l = tempy2l + displ(2,iglob)*hp2
+ tempz2l = tempz2l + displ(3,iglob)*hp2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec)
+ tempx3l = tempx3l + displ(1,iglob)*hp3
+ tempy3l = tempy3l + displ(2,iglob)*hp3
+ tempz3l = tempz3l + displ(3,iglob)*hp3
+ enddo
+
+! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ jacobianl = jacobian(i,j,k,ispec)
+
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ kappal = kappastore(i,j,k,ispec)
+ mul = mustore(i,j,k,ispec)
+
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
+
+! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+
+ sigma_xy = mul*duxdyl_plus_duydxl
+ sigma_xz = mul*duzdxl_plus_duxdzl
+ sigma_yz = mul*duzdyl_plus_duydzl
+
+! form dot product with test vector, symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+
+ enddo
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0.
+ tempy1l = 0.
+ tempz1l = 0.
+
+ tempx2l = 0.
+ tempy2l = 0.
+ tempz2l = 0.
+
+ tempx3l = 0.
+ tempy3l = 0.
+ tempz3l = 0.
+
+ do l=1,NGLLX
+ fac1 = hprimewgll_xx(l,i)
+ tempx1l = tempx1l + tempx1(l,j,k)*fac1
+ tempy1l = tempy1l + tempy1(l,j,k)*fac1
+ tempz1l = tempz1l + tempz1(l,j,k)*fac1
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ fac2 = hprimewgll_yy(l,j)
+ tempx2l = tempx2l + tempx2(i,l,k)*fac2
+ tempy2l = tempy2l + tempy2(i,l,k)*fac2
+ tempz2l = tempz2l + tempz2(i,l,k)*fac2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ fac3 = hprimewgll_zz(l,k)
+ tempx3l = tempx3l + tempx3(i,j,l)*fac3
+ tempy3l = tempy3l + tempy3(i,j,l)*fac3
+ tempz3l = tempz3l + tempz3(i,j,l)*fac3
+ enddo
+
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+
+! sum contributions from each element to the global mesh
+
+ iglob = ibool(i,j,k,ispec)
+
+ accel(1,iglob) = accel(1,iglob) - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
+ accel(2,iglob) = accel(2,iglob) - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
+ accel(3,iglob) = accel(3,iglob) - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
+
+ enddo
+ enddo
+ enddo
+
+ endif ! if (ispec_is_inner(ispec) == phase_is_inner)
+
+ enddo ! spectral element loop
+
+! adding source
+ do isource = 1,NSOURCES
+
+ if (ispec_is_inner(ispec_selected_source(isource)) .eqv. phase_is_inner) then
+
+ if(FASTER_SOURCES_POINTS_ONLY) then
+
+! add the source (only if this proc carries the source)
+ if(myrank == islice_selected_source(isource)) then
+
+ iglob = ibool(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec_selected_source(isource))
+ f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
+ t0 = 1.2d0/f0
+if (it == 1 .and. myrank == 0) then
+
+print *,'using a source of dominant frequency ',f0
+print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+endif
+
+
+ ! we use nu_source(:,3) here because we want a source normal to the surface.
+ ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
+ !accel(:,iglob) = accel(:,iglob) + &
+ ! sngl(nu_source(:,3,isource) * 10000000.d0 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
+ ! exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
+accel(:,iglob) = accel(:,iglob) + &
+ sngl(nu_source(:,3,isource) * 1.d10 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
+ exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
+
+ endif
+ endif
+
+ endif
+
+ enddo
+
+end subroutine compute_forces
Modified: seismo/3D/SPECFEM3D/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/constants.h.in 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/constants.h.in 2008-10-31 19:04:49 UTC (rev 13206)
@@ -116,6 +116,42 @@
character(len=150), parameter :: LOCAL_PATH_Q = '/ibrixfs1/scratch/lqy/DATABASES_MPI_Q/'
!------------------------------------------------------
+! nlegoff -- Variables that should be read/computed elsewhere.
+! Temporarily declared here.
+!------------------------------------------------------
+! whether or not an external mesh is used (provided by CUBIT for example)
+ logical, parameter :: USE_EXTERNAL_MESH = .false.
+
+! no lagrange interpolation on seismograms (we take the value on one NGLL point)
+ logical, parameter :: FASTER_RECEIVERS_POINTS_ONLY = .false.
+ logical, parameter :: FASTER_SOURCES_POINTS_ONLY = .false.
+
+! the receivers can be located inside the model
+ logical, parameter :: RECVS_CAN_BE_BURIED_EXT_MESH = .true.
+ logical, parameter :: SOURCES_CAN_BE_BURIED_EXT_MESH = .true.
+
+! the seismograms are normal to surface
+ logical, parameter :: EXT_MESH_RECV_NORMAL = .false.
+
+!
+ logical, parameter :: EXTERNAL_MESH_MOVIE_SURFACE = .false.
+ logical, parameter :: EXTERNAL_MESH_CREATE_SHAKEMAP = .false.
+
+! deltat
+ double precision, parameter :: DT_ext_mesh = 0.001d0
+
+! NSTEP
+ integer, parameter :: NSTEP_ext_mesh = 100
+
+! number of nodes per element as provided by the external mesh
+ integer, parameter :: ESIZE = 8
+
+! geometry tolerance parameter to calculate number of independent grid points
+! sensitive to actual size of model, assumes reference sphere of radius 1
+! this is an absolute value for normalized coordinates in the Earth
+ double precision, parameter :: SMALLVAL_TOL = 1.d-10
+
+!------------------------------------------------------
!----------- do not modify anything below -------------
!------------------------------------------------------
Modified: seismo/3D/SPECFEM3D/trunk/create_movie_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/create_movie_AVS_DX.f90 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/create_movie_AVS_DX.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -40,13 +40,13 @@
! coefficient of power law used for non linear scaling
logical, parameter :: NONLINEAR_SCALING = .true.
- real(kind=CUSTOM_REAL), parameter :: POWER_SCALING = 0.25_CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: POWER_SCALING = 0.13_CUSTOM_REAL
integer it,it1,it2,ivalue,nspectot_AVS_max,ispec
integer iformat,nframes,iframe,inumber,inorm,iscaling_shake
integer ibool_number,ibool_number1,ibool_number2,ibool_number3,ibool_number4
- logical USE_OPENDX,USE_AVS,USE_GMT,UNIQUE_FILE,plot_shaking_map
+ logical USE_OPENDX,USE_AVS,UNIQUE_FILE,plot_shaking_map
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: x,y,z,display
real(kind=CUSTOM_REAL) xcoord,ycoord,zcoord
@@ -90,7 +90,6 @@
OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL, &
BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS,SAVE_FORWARD
logical ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
- double precision zscaling
character(len=150) OUTPUT_FILES,LOCAL_PATH,MODEL
@@ -104,6 +103,16 @@
NSPEC2D_BOTTOM,NSPEC2D_TOP, &
NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB
+!!!! NL NL for external meshes
+ real(kind=CUSTOM_REAL), parameter :: RADIUS_TO_MUTE = 1000._CUSTOM_REAL
+ logical, parameter :: MUTE_SOURCE = .true.
+ real(kind=CUSTOM_REAL), parameter :: X_SOURCE_EXT_MESH = -9023.021484375
+ real(kind=CUSTOM_REAL), parameter :: Y_SOURCE_EXT_MESH = 6123.611328125
+ real(kind=CUSTOM_REAL), parameter :: Z_SOURCE_EXT_MESH = 17.96331405639648
+ integer, parameter :: NSPEC_SURFACE_EXT_MESH = 15808*4
+
+!!!! NL NL
+
! ************** PROGRAM STARTS HERE **************
print *
@@ -151,10 +160,19 @@
endif
print *
- if(USE_HIGHRES_FOR_MOVIES) then
- ilocnum = NGLLSQUARE*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ if (USE_EXTERNAL_MESH) then
+ NPROC = 1
+ if(USE_HIGHRES_FOR_MOVIES) then
+ ilocnum = NSPEC_SURFACE_EXT_MESH*NGLLSQUARE
+ else
+ ilocnum = NSPEC_SURFACE_EXT_MESH*NGNOD2D_AVS_DX
+ endif
else
- ilocnum = NGNOD2D_AVS_DX*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ if(USE_HIGHRES_FOR_MOVIES) then
+ ilocnum = NGLLSQUARE*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ else
+ ilocnum = NGNOD2D_AVS_DX*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ endif
endif
allocate(store_val_x(ilocnum,0:NPROC-1))
allocate(store_val_y(ilocnum,0:NPROC-1))
@@ -162,6 +180,7 @@
allocate(store_val_ux(ilocnum,0:NPROC-1))
allocate(store_val_uy(ilocnum,0:NPROC-1))
allocate(store_val_uz(ilocnum,0:NPROC-1))
+
print *,'1 = create files in OpenDX format'
print *,'2 = create files in AVS UCD format with individual files'
@@ -175,33 +194,21 @@
if(iformat == 1) then
USE_OPENDX = .true.
USE_AVS = .false.
- USE_GMT = .false.
UNIQUE_FILE = .false.
else if(iformat == 2) then
USE_OPENDX = .false.
USE_AVS = .true.
- USE_GMT = .false.
UNIQUE_FILE = .false.
else if(iformat == 3) then
USE_OPENDX = .false.
USE_AVS = .true.
- USE_GMT = .false.
UNIQUE_FILE = .true.
else
USE_OPENDX = .false.
USE_AVS = .false.
- USE_GMT = .true.
UNIQUE_FILE = .false.
endif
- if(.not. USE_GMT) then
- print *
- print *,'scaling to apply to Z to amplify topography (1. to do nothing, 0. to get flat surface):'
- read(5,*) zscaling
- else
- zscaling = 0.
- endif
-
print *,'movie frames have been saved every ',NTSTEP_BETWEEN_FRAMES,' time steps'
print *
@@ -262,10 +269,18 @@
endif
! define the total number of elements at the surface
- if(USE_HIGHRES_FOR_MOVIES) then
- nspectot_AVS_max = NEX_XI * NEX_ETA * (NGLLX-1) * (NGLLY-1)
+ if (USE_EXTERNAL_MESH) then
+ if(USE_HIGHRES_FOR_MOVIES) then
+ nspectot_AVS_max = NSPEC_SURFACE_EXT_MESH * (NGLLX-1) * (NGLLY-1)
+ else
+ nspectot_AVS_max = NSPEC_SURFACE_EXT_MESH
+ endif
else
- nspectot_AVS_max = NEX_XI * NEX_ETA
+ if(USE_HIGHRES_FOR_MOVIES) then
+ nspectot_AVS_max = NEX_XI * NEX_ETA * (NGLLX-1) * (NGLLY-1)
+ else
+ nspectot_AVS_max = NEX_XI * NEX_ETA
+ endif
endif
print *
@@ -340,7 +355,8 @@
! reset point number
ipoin = 0
- do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ !do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ do ispecloc = 1, NSPEC_SURFACE_EXT_MESH
if(USE_HIGHRES_FOR_MOVIES) then
! assign the OpenDX "elements"
@@ -354,9 +370,6 @@
ycoord = store_val_y(ipoin,iproc)
zcoord = store_val_z(ipoin,iproc)
-! amplify topography, otherwise too flat to see anything
- zcoord = zcoord * zscaling
-
vectorx = store_val_ux(ipoin,iproc)
vectory = store_val_uy(ipoin,iproc)
vectorz = store_val_uz(ipoin,iproc)
@@ -366,6 +379,15 @@
z(i,j) = zcoord
if(plot_shaking_map) then
+!!!! NL NL mute value near source
+ if ( (sqrt(((x(i,j) - (X_SOURCE_EXT_MESH))**2 + &
+ (y(i,j) - (Y_SOURCE_EXT_MESH))**2 + &
+ (z(i,j) - (Z_SOURCE_EXT_MESH))**2)) < RADIUS_TO_MUTE) &
+ .and. MUTE_SOURCE) then
+
+ display(i,j) = 0.
+ else
+
if(inorm == 1) then
display(i,j) = vectorx
else if(inorm == 2) then
@@ -373,8 +395,13 @@
else
display(i,j) = vectorz
endif
+ endif
else
- display(i,j) = vectorz
+ if (USE_EXTERNAL_MESH) then
+ display(i,j) = sqrt(vectorz**2+vectory**2+vectorx**2)
+ else
+ display(i,j) = vectorz
+ endif
endif
enddo
@@ -429,9 +456,6 @@
ycoord = store_val_y(ipoin,iproc)
zcoord = store_val_z(ipoin,iproc)
-! amplify topography, otherwise too flat to see anything
- zcoord = zcoord * zscaling
-
vectorx = store_val_ux(ipoin,iproc)
vectory = store_val_uy(ipoin,iproc)
vectorz = store_val_uz(ipoin,iproc)
@@ -444,6 +468,16 @@
! or show norm of vector if shaking map
! for shaking map, norm of U stored in ux, V in uy and A in uz
if(plot_shaking_map) then
+!!!! NL NL mute value near source
+ if ( (sqrt(((x(i,j) - (X_SOURCE_EXT_MESH))**2 + &
+ (y(i,j) - (Y_SOURCE_EXT_MESH))**2 + &
+ (z(i,j) - (Z_SOURCE_EXT_MESH))**2)) < RADIUS_TO_MUTE) &
+ .and. MUTE_SOURCE) then
+
+ display(i,j) = 0.
+ else
+
+
if(inorm == 1) then
field_display(ilocnum+ieoff) = dble(vectorx)
else if(inorm == 2) then
@@ -451,8 +485,13 @@
else
field_display(ilocnum+ieoff) = dble(vectorz)
endif
+ endif
else
- field_display(ilocnum+ieoff) = dble(vectorz)
+ if (USE_EXTERNAL_MESH) then
+ field_display(ilocnum+ieoff) =sqrt(vectorz**2+vectory**2+vectorx**2)
+ else
+ field_display(ilocnum+ieoff) = dble(vectorz)
+ endif
endif
enddo
@@ -469,7 +508,12 @@
!--- sort the list based upon coordinates to get rid of multiples
print *,'sorting list of points'
- call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot,UTM_X_MIN,UTM_X_MAX)
+ if (USE_EXTERNAL_MESH) then
+ call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot, &
+ dble(minval(store_val_x(:,0))),dble(maxval(store_val_x(:,0))))
+ else
+ call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot,UTM_X_MIN,UTM_X_MAX)
+ endif
!--- print total number of points found
print *
@@ -501,7 +545,6 @@
print *
print *,'minimum amplitude in current snapshot after removal = ',min_field_current
print *,'maximum amplitude in current snapshot after removal = ',max_field_current
- if(inorm == 3) print *,'maximum corresponds to ',max_field_current/9.81,' g'
print *
endif
@@ -581,9 +624,6 @@
write(outputname,"('/AVS_shaking_map.inp')")
open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
- else if(USE_GMT) then
- write(outputname,"('/gmt_shaking_map.xyz')")
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
else
stop 'wrong output format selected'
endif
@@ -606,33 +646,14 @@
open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
endif
- else if(USE_GMT) then
- write(outputname,"('/gmt_movie_',i6.6,'.xyz')") ivalue
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
else
stop 'wrong output format selected'
endif
endif
- if(USE_GMT) then
+ if(.false.) then
-! output list of points
- mask_point = .false.
- do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
-! four points for each element
- do ilocnum = 1,NGNOD2D_AVS_DX
- ibool_number = iglob(ilocnum+ieoff)
- if(.not. mask_point(ibool_number)) then
- call utm_geo(long,lat,xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff), &
- UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
- write(11,*) long,lat,field_display(ilocnum+ieoff)
- endif
- mask_point(ibool_number) = .true.
- enddo
- enddo
-
else
! if unique file, output geometry only once
@@ -734,7 +755,7 @@
endif
mask_point(ibool_number) = .true.
enddo
- enddo
+ enddo
! define OpenDX field
if(USE_OPENDX) then
@@ -762,7 +783,7 @@
print *
if(USE_OPENDX) print *,'DX files are stored in ', trim(OUTPUT_FILES), '/DX_*.dx'
if(USE_AVS) print *,'AVS files are stored in ', trim(OUTPUT_FILES), '/AVS_*.inp'
- if(USE_GMT) print *,'GMT files are stored in ', trim(OUTPUT_FILES), '/gmt_*.xyz'
+
print *
@@ -825,7 +846,10 @@
double precision UTM_X_MIN,UTM_X_MAX
! define geometrical tolerance based upon typical size of the model
- SMALLVALTOL = 1.d-10 * dabs(UTM_X_MAX - UTM_X_MIN)
+ SMALLVALTOL = 1.d-10 * dabs(UTM_X_MAX - UTM_X_MIN)
+ print *, 'UTM_X_MAX', UTM_X_MAX
+ print *, 'UTM_X_MIN', UTM_X_MIN
+ print *, 'SMALLVALTOL', SMALLVALTOL
! dynamically allocate arrays
allocate(ind(npointot))
Modified: seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -1153,3 +1153,972 @@
end subroutine create_regions_mesh
+!
+!----
+!
+
+ subroutine create_regions_mesh_ext_mesh(ibool, &
+ xstore,ystore,zstore,npx,npy,iproc_xi,iproc_eta,nspec, &
+ volume_local,area_local_bottom,area_local_top, &
+ NGLOB_AB,npointot, &
+ myrank,LOCAL_PATH, &
+ nnodes_ext_mesh,nelmnts_ext_mesh, &
+ nodes_coords_ext_mesh,elmnts_ext_mesh,mat_ext_mesh, &
+ ninterface_ext_mesh,max_interface_size_ext_mesh, &
+ my_neighbours_ext_mesh,my_nelmnts_neighbours_ext_mesh,my_interfaces_ext_mesh, &
+ ibool_interfaces_ext_mesh,nibool_interfaces_ext_mesh &
+ )
+
+! create the different regions of the mesh
+
+ implicit none
+
+ include "constants.h"
+
+! number of spectral elements in each block
+ integer nspec
+
+ integer npx,npy
+ integer npointot
+
+ character(len=150) LOCAL_PATH
+
+! arrays with the mesh
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
+
+ double precision xstore_local(NGLLX,NGLLY,NGLLZ)
+ double precision ystore_local(NGLLX,NGLLY,NGLLZ)
+ double precision zstore_local(NGLLX,NGLLY,NGLLZ)
+
+ double precision xmesh,ymesh,zmesh
+
+ integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+! data from the external mesh
+ integer :: nnodes_ext_mesh,nelmnts_ext_mesh
+ double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
+ integer, dimension(ESIZE,nelmnts_ext_mesh) :: elmnts_ext_mesh
+ integer, dimension(nelmnts_ext_mesh) :: mat_ext_mesh
+ double precision, external :: materials_ext_mesh
+ integer :: ninterface_ext_mesh,max_interface_size_ext_mesh
+ integer, dimension(ninterface_ext_mesh) :: my_neighbours_ext_mesh
+ integer, dimension(ninterface_ext_mesh) :: my_nelmnts_neighbours_ext_mesh
+ integer, dimension(6,max_interface_size_ext_mesh,ninterface_ext_mesh) :: my_interfaces_ext_mesh
+ integer, dimension(NGLLX*NGLLX*max_interface_size_ext_mesh,ninterface_ext_mesh) :: ibool_interfaces_ext_mesh
+ integer, dimension(ninterface_ext_mesh) :: nibool_interfaces_ext_mesh
+
+! for MPI buffers
+ integer, dimension(:), allocatable :: reorder_interface_ext_mesh,ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh
+ integer, dimension(:), allocatable :: nibool_interfaces_ext_mesh_true
+ integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh_dummy
+ integer, dimension(:), allocatable :: ibool_interface_ext_mesh_dummy
+ double precision, dimension(:), allocatable :: work_ext_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: xstore_dummy
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: ystore_dummy
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: zstore_dummy
+
+! Gauss-Lobatto-Legendre points and weights of integration
+ double precision, dimension(:), allocatable :: xigll,yigll,zigll,wxgll,wygll,wzgll
+
+! 3D shape functions and their derivatives
+ double precision, dimension(:,:,:,:), allocatable :: shape3D
+ double precision, dimension(:,:,:,:,:), allocatable :: dershape3D
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+
+! the jacobian
+ real(kind=CUSTOM_REAL) jacobianl
+
+! arrays with mesh parameters
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: xixstore,xiystore,xizstore, &
+ etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,jacobianstore
+
+! for model density
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore
+
+! proc numbers for MPI
+ integer myrank
+
+! check area and volume of the final mesh
+ double precision weight
+ double precision area_local_bottom,area_local_top
+ double precision volume_local
+
+! variables for creating array ibool (some arrays also used for AVS or DX files)
+ integer, dimension(:), allocatable :: iglob,locval
+ logical, dimension(:), allocatable :: ifseg
+ double precision, dimension(:), allocatable :: xp,yp,zp
+
+ integer nglob,NGLOB_AB
+ integer ieoff,ilocnum
+ integer ier
+ integer iinterface
+
+! mass matrix
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass
+
+! ---------------------------
+
+! name of the database file
+ character(len=150) prname
+
+ integer i,j,k,ia,ispec,iglobnum,itype_element
+ integer iproc_xi,iproc_eta
+
+ double precision rho,vp,vs
+ double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+
+! for the harvard 3D salton sea model
+ double precision :: umesh, vmesh, wmesh, vp_st, vs_st, rho_st
+
+! mask to sort ibool
+ integer, dimension(:), allocatable :: mask_ibool
+ integer, dimension(:,:,:,:), allocatable :: copy_ibool_ori
+ integer :: inumber
+
+! **************
+
+! create the name for the database of the current slide and region
+ call create_name_database(prname,myrank,LOCAL_PATH)
+
+! Gauss-Lobatto-Legendre points of integration
+ allocate(xigll(NGLLX))
+ allocate(yigll(NGLLY))
+ allocate(zigll(NGLLZ))
+
+! Gauss-Lobatto-Legendre weights of integration
+ allocate(wxgll(NGLLX))
+ allocate(wygll(NGLLY))
+ allocate(wzgll(NGLLZ))
+
+! 3D shape functions and their derivatives
+ allocate(shape3D(NGNOD,NGLLX,NGLLY,NGLLZ))
+ allocate(dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ))
+
+! array with model density
+ allocate(rhostore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(kappastore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(mustore(NGLLX,NGLLY,NGLLZ,nspec))
+
+! arrays with mesh parameters
+ allocate(xixstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(xiystore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(xizstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(etaxstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(etaystore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(etazstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(gammaxstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(gammaystore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(gammazstore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(jacobianstore(NGLLX,NGLLY,NGLLZ,nspec))
+
+! set up coordinates of the Gauss-Lobatto-Legendre points
+ call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
+
+! if number of points is odd, the middle abscissa is exactly zero
+ if(mod(NGLLX,2) /= 0) xigll((NGLLX-1)/2+1) = ZERO
+ if(mod(NGLLY,2) /= 0) yigll((NGLLY-1)/2+1) = ZERO
+ if(mod(NGLLZ,2) /= 0) zigll((NGLLZ-1)/2+1) = ZERO
+
+! get the 3-D shape functions
+ call get_shape3D(myrank,shape3D,dershape3D,xigll,yigll,zigll)
+
+! allocate memory for arrays
+ allocate(iglob(npointot))
+ allocate(locval(npointot))
+ allocate(ifseg(npointot))
+ allocate(xp(npointot))
+ allocate(yp(npointot))
+ allocate(zp(npointot))
+
+!---
+
+ xstore(:,:,:,:) = 0.d0
+ ystore(:,:,:,:) = 0.d0
+ zstore(:,:,:,:) = 0.d0
+
+ do ispec = 1, nspec
+ !call get_xyzelm(xelm, yelm, zelm, ispec, elmnts_ext_mesh, nodes_coords_ext_mesh, nspec, nnodes_ext_mesh)
+ do ia = 1,NGNOD
+ xelm(ia) = nodes_coords_ext_mesh(1,elmnts_ext_mesh(ia,ispec))
+ yelm(ia) = nodes_coords_ext_mesh(2,elmnts_ext_mesh(ia,ispec))
+ zelm(ia) = nodes_coords_ext_mesh(3,elmnts_ext_mesh(ia,ispec))
+ enddo
+
+ call calc_jacobian(myrank,xixstore,xiystore,xizstore, &
+ etaxstore,etaystore,etazstore, &
+ gammaxstore,gammaystore,gammazstore,jacobianstore, &
+ xstore,ystore,zstore, &
+ xelm,yelm,zelm,shape3D,dershape3D,ispec,nspec)
+
+ enddo
+
+! kappastore and mustore
+ do ispec = 1, nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ kappastore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
+ (materials_ext_mesh(2,mat_ext_mesh(ispec))*materials_ext_mesh(2,mat_ext_mesh(ispec)) - &
+ 4.d0*materials_ext_mesh(3,mat_ext_mesh(ispec))*materials_ext_mesh(3,mat_ext_mesh(ispec))/3.d0)
+ mustore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))*materials_ext_mesh(3,mat_ext_mesh(ispec))*&
+ materials_ext_mesh(3,mat_ext_mesh(ispec))
+ enddo
+ enddo
+ enddo
+ enddo
+
+ locval = 0
+ ifseg = .false.
+ xp = 0.d0
+ yp = 0.d0
+ zp = 0.d0
+
+ do ispec=1,nspec
+ ieoff = NGLLX * NGLLY * NGLLZ * (ispec-1)
+ ilocnum = 0
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ilocnum = ilocnum + 1
+ xp(ilocnum+ieoff) = xstore(i,j,k,ispec)
+ yp(ilocnum+ieoff) = ystore(i,j,k,ispec)
+ zp(ilocnum+ieoff) = zstore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ call get_global(nspec,xp,yp,zp,ibool,locval,ifseg,nglob,npointot, &
+ minval(nodes_coords_ext_mesh(1,:)),maxval(nodes_coords_ext_mesh(1,:)))
+
+ deallocate(xp,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(yp,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(zp,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(locval,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(ifseg,stat=ier); if(ier /= 0) stop 'error in deallocate'
+
+!
+!- we can create a new indirect addressing to reduce cache misses
+!
+ allocate(copy_ibool_ori(NGLLX,NGLLY,NGLLZ,nspec),stat=ier); if(ier /= 0) stop 'error in allocate'
+ allocate(mask_ibool(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
+
+ mask_ibool(:) = -1
+ copy_ibool_ori(:,:,:,:) = ibool(:,:,:,:)
+
+ inumber = 0
+ do ispec=1,nspec
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ if(mask_ibool(copy_ibool_ori(i,j,k,ispec)) == -1) then
+! create a new point
+ inumber = inumber + 1
+ ibool(i,j,k,ispec) = inumber
+ mask_ibool(copy_ibool_ori(i,j,k,ispec)) = inumber
+ else
+! use an existing point created previously
+ ibool(i,j,k,ispec) = mask_ibool(copy_ibool_ori(i,j,k,ispec))
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+
+ deallocate(copy_ibool_ori,stat=ier); if(ier /= 0) stop 'error in deallocate'
+ deallocate(mask_ibool,stat=ier); if(ier /= 0) stop 'error in deallocate'
+
+ allocate(xstore_dummy(nglob))
+ allocate(ystore_dummy(nglob))
+ allocate(zstore_dummy(nglob))
+ do ispec = 1, nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglobnum = ibool(i,j,k,ispec)
+ xstore_dummy(iglobnum) = xstore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ do ispec = 1, nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglobnum = ibool(i,j,k,ispec)
+ ystore_dummy(iglobnum) = ystore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ do ispec = 1, nspec
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglobnum = ibool(i,j,k,ispec)
+ zstore_dummy(iglobnum) = zstore(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+
+! creating mass matrix (will be fully assembled with MPI in the solver)
+ allocate(rmass(nglob))
+ rmass(:) = 0._CUSTOM_REAL
+
+ do ispec=1,nspec
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ weight=wxgll(i)*wygll(j)*wzgll(k)
+ iglobnum=ibool(i,j,k,ispec)
+
+ jacobianl=jacobianstore(i,j,k,ispec)
+
+! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass(iglobnum) = rmass(iglobnum) + &
+ sngl(dble(materials_ext_mesh(1,mat_ext_mesh(ispec))) * dble(jacobianl) * weight)
+ else
+ rmass(iglobnum) = rmass(iglobnum) + materials_ext_mesh(1,mat_ext_mesh(ispec)) * jacobianl * weight
+ endif
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ call prepare_assemble_MPI (nelmnts_ext_mesh,ibool, &
+ elmnts_ext_mesh, ESIZE, &
+ nglob, &
+ ninterface_ext_mesh, max_interface_size_ext_mesh, &
+ my_nelmnts_neighbours_ext_mesh, my_interfaces_ext_mesh, &
+ ibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh &
+ )
+
+! sort ibool comm buffers lexicographically
+ allocate(nibool_interfaces_ext_mesh_true(ninterface_ext_mesh))
+
+ do iinterface = 1, ninterface_ext_mesh
+
+ allocate(xp(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(yp(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(zp(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(locval(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ifseg(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(reorder_interface_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ibool_interface_ext_mesh_dummy(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ind_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(ninseg_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(iwork_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(work_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+
+ do ilocnum = 1, nibool_interfaces_ext_mesh(iinterface)
+ xp(ilocnum) = xstore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
+ yp(ilocnum) = ystore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
+ zp(ilocnum) = zstore_dummy(ibool_interfaces_ext_mesh(ilocnum,iinterface))
+ enddo
+
+ call sort_array_coordinates(nibool_interfaces_ext_mesh(iinterface),xp,yp,zp, &
+ ibool_interfaces_ext_mesh(1:nibool_interfaces_ext_mesh(iinterface),iinterface), &
+ reorder_interface_ext_mesh,locval,ifseg,nibool_interfaces_ext_mesh_true(iinterface), &
+ ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh,work_ext_mesh)
+
+ deallocate(xp)
+ deallocate(yp)
+ deallocate(zp)
+ deallocate(locval)
+ deallocate(ifseg)
+ deallocate(reorder_interface_ext_mesh)
+ deallocate(ibool_interface_ext_mesh_dummy)
+ deallocate(ind_ext_mesh)
+ deallocate(ninseg_ext_mesh)
+ deallocate(iwork_ext_mesh)
+ deallocate(work_ext_mesh)
+
+ enddo
+
+! save the binary files
+ call create_name_database(prname,myrank,LOCAL_PATH)
+ open(unit=IOUT,file=prname(1:len_trim(prname))//'external_mesh.bin',status='unknown',action='write',form='unformatted')
+ write(IOUT) nspec
+ write(IOUT) nglob
+
+ write(IOUT) xixstore
+ write(IOUT) xiystore
+ write(IOUT) xizstore
+ write(IOUT) etaxstore
+ write(IOUT) etaystore
+ write(IOUT) etazstore
+ write(IOUT) gammaxstore
+ write(IOUT) gammaystore
+ write(IOUT) gammazstore
+
+ write(IOUT) jacobianstore
+ write(IOUT) kappastore
+ write(IOUT) mustore
+
+ write(IOUT) rmass
+
+ write(IOUT) ibool
+
+ write(IOUT) xstore_dummy
+ write(IOUT) ystore_dummy
+ write(IOUT) zstore_dummy
+
+ write(IOUT) ninterface_ext_mesh
+ write(IOUT) maxval(nibool_interfaces_ext_mesh)
+ write(IOUT) my_neighbours_ext_mesh
+ write(IOUT) nibool_interfaces_ext_mesh
+ allocate(ibool_interfaces_ext_mesh_dummy(maxval(nibool_interfaces_ext_mesh),ninterface_ext_mesh))
+ do i = 1, ninterface_ext_mesh
+ ibool_interfaces_ext_mesh_dummy = ibool_interfaces_ext_mesh(1:maxval(nibool_interfaces_ext_mesh),:)
+ enddo
+ write(IOUT) ibool_interfaces_ext_mesh_dummy
+ close(IOUT)
+
+ end subroutine create_regions_mesh_ext_mesh
+
+!
+!----
+!
+
+ double precision function materials_ext_mesh(i,j)
+
+ implicit none
+
+ integer :: i,j
+
+ select case (j)
+ case (1)
+ select case (i)
+ case (1)
+ materials_ext_mesh = 2000.d0
+ case (2)
+ materials_ext_mesh = 3000.d0
+ case (3)
+ materials_ext_mesh = 1732.051d0
+ case default
+ call stop_all()
+ end select
+ case (2)
+ select case (i)
+ case (1)
+ materials_ext_mesh = 2000.d0
+ case (2)
+ materials_ext_mesh = 900.d0
+ case (3)
+ materials_ext_mesh = 500.d0
+ case default
+ call stop_all()
+ end select
+ case default
+ call stop_all()
+ end select
+
+ end function materials_ext_mesh
+
+!
+!----
+!
+
+subroutine prepare_assemble_MPI (nelmnts,ibool, &
+ knods, ngnode, &
+ npoin, &
+ ninterface, max_interface_size, &
+ my_nelmnts_neighbours, my_interfaces, &
+ ibool_interfaces_asteroid, &
+ nibool_interfaces_asteroid &
+ )
+
+ implicit none
+
+ include 'constants.h'
+
+ integer, intent(in) :: nelmnts, npoin, ngnode
+ integer, dimension(ngnode,nelmnts), intent(in) :: knods
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nelmnts), intent(in) :: ibool
+
+ integer :: ninterface
+ integer :: max_interface_size
+ integer, dimension(ninterface) :: my_nelmnts_neighbours
+ integer, dimension(6,max_interface_size,ninterface) :: my_interfaces
+ integer, dimension(NGLLX*NGLLX*max_interface_size,ninterface) :: &
+ ibool_interfaces_asteroid
+ integer, dimension(ninterface) :: &
+ nibool_interfaces_asteroid
+
+ integer :: num_interface
+ integer :: ispec_interface
+
+ logical, dimension(npoin) :: mask_ibool_asteroid
+
+ integer :: ixmin, ixmax
+ integer :: iymin, iymax
+ integer :: izmin, izmax
+ integer, dimension(ngnode) :: n
+ integer :: e1, e2, e3, e4
+ integer :: type
+ integer :: ispec
+
+ integer :: k
+ integer :: npoin_interface_asteroid
+
+ integer :: ix,iy,iz
+
+
+ ibool_interfaces_asteroid(:,:) = 0
+ nibool_interfaces_asteroid(:) = 0
+
+ do num_interface = 1, ninterface
+ npoin_interface_asteroid = 0
+ mask_ibool_asteroid(:) = .false.
+
+ do ispec_interface = 1, my_nelmnts_neighbours(num_interface)
+ ispec = my_interfaces(1,ispec_interface,num_interface)
+ type = my_interfaces(2,ispec_interface,num_interface)
+ do k = 1, ngnode
+ n(k) = knods(k,ispec)
+ end do
+ e1 = my_interfaces(3,ispec_interface,num_interface)
+ e2 = my_interfaces(4,ispec_interface,num_interface)
+ e3 = my_interfaces(5,ispec_interface,num_interface)
+ e4 = my_interfaces(6,ispec_interface,num_interface)
+ call get_edge(ngnode, n, type, e1, e2, e3, e4, ixmin, ixmax, iymin, iymax, izmin, izmax)
+
+ do iz = min(izmin,izmax), max(izmin,izmax)
+ do iy = min(iymin,iymax), max(iymin,iymax)
+ do ix = min(ixmin,ixmax), max(ixmin,ixmax)
+
+ if(.not. mask_ibool_asteroid(ibool(ix,iy,iz,ispec))) then
+ mask_ibool_asteroid(ibool(ix,iy,iz,ispec)) = .true.
+ npoin_interface_asteroid = npoin_interface_asteroid + 1
+ ibool_interfaces_asteroid(npoin_interface_asteroid,num_interface)=&
+ ibool(ix,iy,iz,ispec)
+ end if
+ end do
+ end do
+ end do
+
+ end do
+ nibool_interfaces_asteroid(num_interface) = npoin_interface_asteroid
+
+
+ end do
+
+end subroutine prepare_assemble_MPI
+
+!
+!----
+!
+
+subroutine get_edge ( ngnode, n, type, e1, e2, e3, e4, ixmin, ixmax, iymin, iymax, izmin, izmax )
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: ngnode
+ integer, dimension(ngnode), intent(in) :: n
+ integer, intent(in) :: type, e1, e2, e3, e4
+ integer, intent(out) :: ixmin, ixmax, iymin, iymax, izmin, izmax
+
+ integer, dimension(4) :: en
+ integer :: valence, i
+
+ if ( type == 1 ) then
+ if ( e1 == n(1) ) then
+ ixmin = 1
+ ixmax = 1
+ iymin = 1
+ iymax = 1
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(2) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ iymin = 1
+ iymax = 1
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(3) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ iymin = NGLLY
+ iymax = NGLLY
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(4) ) then
+ ixmin = 1
+ ixmax = 1
+ iymin = NGLLY
+ iymax = NGLLY
+ izmin = 1
+ izmax = 1
+ end if
+ if ( e1 == n(5) ) then
+ ixmin = 1
+ ixmax = 1
+ iymin = 1
+ iymax = 1
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ if ( e1 == n(6) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ iymin = 1
+ iymax = 1
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ if ( e1 == n(7) ) then
+ ixmin = NGLLX
+ ixmax = NGLLX
+ iymin = NGLLY
+ iymax = NGLLY
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ if ( e1 == n(8) ) then
+ ixmin = 1
+ ixmax = 1
+ iymin = NGLLY
+ iymax = NGLLY
+ izmin = NGLLZ
+ izmax = NGLLZ
+ end if
+ else
+ if ( type == 2 ) then
+ if ( e1 == n(1) ) then
+ ixmin = 1
+ iymin = 1
+ izmin = 1
+ if ( e2 == n(2) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(4) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(5) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(2) ) then
+ ixmin = NGLLX
+ iymin = 1
+ izmin = 1
+ if ( e2 == n(3) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(1) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(6) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = NGLLZ
+ end if
+
+ end if
+ if ( e1 == n(3) ) then
+ ixmin = NGLLX
+ iymin = NGLLY
+ izmin = 1
+ if ( e2 == n(4) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(2) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(7) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(4) ) then
+ ixmin = 1
+ iymin = NGLLY
+ izmin = 1
+ if ( e2 == n(1) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(3) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(8) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(5) ) then
+ ixmin = 1
+ iymin = 1
+ izmin = NGLLZ
+ if ( e2 == n(1) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(6) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = NGLLZ
+ end if
+ if ( e2 == n(8) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(6) ) then
+ ixmin = NGLLX
+ iymin = 1
+ izmin = NGLLZ
+ if ( e2 == n(2) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = 1
+ end if
+ if ( e2 == n(7) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ if ( e2 == n(5) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(7) ) then
+ ixmin = NGLLX
+ iymin = NGLLY
+ izmin = NGLLZ
+ if ( e2 == n(3) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(8) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ if ( e2 == n(6) ) then
+ ixmax = NGLLX
+ iymax = 1
+ izmax = NGLLZ
+ end if
+ end if
+ if ( e1 == n(8) ) then
+ ixmin = 1
+ iymin = NGLLY
+ izmin = NGLLZ
+ if ( e2 == n(4) ) then
+ ixmax = 1
+ iymax = NGLLY
+ izmax = 1
+ end if
+ if ( e2 == n(5) ) then
+ ixmax = 1
+ iymax = 1
+ izmax = NGLLZ
+ end if
+ if ( e2 == n(7) ) then
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = NGLLZ
+ end if
+ end if
+
+ else
+ if (type == 4) then
+ en(1) = e1
+ en(2) = e2
+ en(3) = e3
+ en(4) = e4
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(1)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(2)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(3)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(4)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = 1
+ iymin = 1
+ izmin = 1
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = 1
+ endif
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(1)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(2)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(5)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(6)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = 1
+ iymin = 1
+ izmin = 1
+ ixmax = NGLLX
+ iymax = 1
+ izmax = NGLLZ
+ endif
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(2)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(3)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(6)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(7)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = NGLLX
+ iymin = 1
+ izmin = 1
+ ixmax = NGLLX
+ iymax = NGLLZ
+ izmax = NGLLZ
+ endif
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(3)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(4)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(7)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(8)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = 1
+ iymin = NGLLY
+ izmin = 1
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = NGLLZ
+ endif
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(1)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(4)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(5)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(8)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = 1
+ iymin = 1
+ izmin = 1
+ ixmax = 1
+ iymax = NGLLY
+ izmax = NGLLZ
+ endif
+
+ valence = 0
+ do i = 1, 4
+ if ( en(i) == n(5)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(6)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(7)) then
+ valence = valence+1
+ endif
+ if ( en(i) == n(8)) then
+ valence = valence+1
+ endif
+ enddo
+ if ( valence == 4 ) then
+ ixmin = 1
+ iymin = 1
+ izmin = NGLLZ
+ ixmax = NGLLX
+ iymax = NGLLY
+ izmax = NGLLZ
+ endif
+
+ else
+ stop 'ERROR get_edge'
+ endif
+
+ end if
+ end if
+
+end subroutine get_edge
Modified: seismo/3D/SPECFEM3D/trunk/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/locate_receivers.f90 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/locate_receivers.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -32,7 +32,9 @@
xi_receiver,eta_receiver,gamma_receiver,station_name,network_name,nu, &
NPROC,utm_x_source,utm_y_source, &
TOPOGRAPHY,itopo_bathy,UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION, &
- NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO)
+ NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO, &
+ iglob_is_surface_external_mesh,ispec_is_surface_external_mesh &
+ )
implicit none
@@ -57,6 +59,12 @@
integer itopo_bathy(NX_TOPO,NY_TOPO)
double precision long_corner,lat_corner,ratio_xi,ratio_eta
+! for surface locating and normal computing with external mesh
+ integer :: pt0_ix,pt0_iy,pt0_iz,pt1_ix,pt1_iy,pt1_iz,pt2_ix,pt2_iy,pt2_iz
+ real(kind=CUSTOM_REAL), dimension(3) :: u_vector,v_vector,w_vector
+ logical, dimension(NGLOB_AB) :: iglob_is_surface_external_mesh
+ logical, dimension(NSPEC_AB) :: ispec_is_surface_external_mesh
+
integer, allocatable, dimension(:) :: ix_initial_guess,iy_initial_guess,iz_initial_guess
integer iprocloop
@@ -68,7 +76,7 @@
double precision, allocatable, dimension(:,:) :: x_found_all,y_found_all,z_found_all
integer irec
- integer i,j,k,ispec,iglob
+ integer i,j,k,ispec,iglob,imin,imax,jmin,jmax,kmin,kmax
integer icornerlong,icornerlat
double precision utm_x_source,utm_y_source
@@ -120,7 +128,8 @@
integer, allocatable, dimension(:,:) :: ispec_selected_rec_all
double precision, allocatable, dimension(:) :: stlat,stlon,stele,stbur,stutm_x,stutm_y
double precision, allocatable, dimension(:,:) :: xi_receiver_all,eta_receiver_all,gamma_receiver_all
-
+ double precision, allocatable, dimension(:,:,:,:) :: nu_all
+
character(len=150) OUTPUT_FILES
! **************
@@ -179,6 +188,7 @@
allocate(y_found_all(nrec,0:NPROC-1))
allocate(z_found_all(nrec,0:NPROC-1))
allocate(final_distance_all(nrec,0:NPROC-1))
+ allocate(nu_all(3,3,nrec,0:NPROC-1))
! loop on all the stations
do irec=1,nrec
@@ -190,7 +200,8 @@
if (ios /= 0) call exit_mpi(myrank, 'Error reading station file '//trim(rec_filename))
! convert station location to UTM
- call utm_geo(stlon(irec),stlat(irec),stutm_x(irec),stutm_y(irec),UTM_PROJECTION_ZONE,ILONGLAT2UTM,SUPPRESS_UTM_PROJECTION)
+ call utm_geo(stlon(irec),stlat(irec),stutm_x(irec),stutm_y(irec),UTM_PROJECTION_ZONE,ILONGLAT2UTM, &
+ (SUPPRESS_UTM_PROJECTION .or. USE_EXTERNAL_MESH))
! compute horizontal distance between source and receiver in km
horiz_dist(irec) = dsqrt((stutm_y(irec)-utm_y_source)**2 + (stutm_x(irec)-utm_x_source)**2) / 1000.
@@ -220,6 +231,7 @@
nu(3,2,irec) = 0.d0
nu(3,3,irec) = 1.d0
+ if (.not. USE_EXTERNAL_MESH) then
! compute elevation of topography at the receiver location
! we assume that receivers are always at the surface i.e. not buried
if(TOPOGRAPHY) then
@@ -265,21 +277,58 @@
z_target(irec) = elevation(irec) - stbur(irec)
if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
+ else
+
+ x_target(irec) = stutm_x(irec)
+ y_target(irec) = stutm_y(irec)
+ z_target(irec) = stbur(irec)
+ if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
+
+ endif ! of if (.not. USE_EXTERNAL_MESH)
+
! examine top of the elements only (receivers always at the surface)
! k = NGLLZ
+ ispec_selected_rec(irec) = 0
+
do ispec=1,NSPEC_AB
-! modification by Qinya Liu: idoubling is no longer used because receivers can now be in depth
-! (in previous versions, receivers were always assumed to be at the surface)
+! define the interval in which we look for points
+ if(FASTER_RECEIVERS_POINTS_ONLY) then
+ imin = 1
+ imax = NGLLX
+ jmin = 1
+ jmax = NGLLY
+
+ kmin = 1
+ kmax = NGLLZ
+
+ else
! loop only on points inside the element
! exclude edges to ensure this point is not shared with other elements
- do k=2,NGLLZ-1
- do j=2,NGLLY-1
- do i=2,NGLLX-1
+ imin = 2
+ imax = NGLLX - 1
+ jmin = 2
+ jmax = NGLLY - 1
+
+ kmin = 2
+ kmax = NGLLZ - 1
+ endif
+
+ do k = kmin,kmax
+ do j = jmin,jmax
+ do i = imin,imax
+
iglob = ibool(i,j,k,ispec)
+
+ if (USE_EXTERNAL_MESH .and. (.not. RECVS_CAN_BE_BURIED_EXT_MESH)) then
+ if ((.not. iglob_is_surface_external_mesh(iglob)) .or. (.not. ispec_is_surface_external_mesh(ispec))) then
+ cycle
+ endif
+ endif
+
dist = dsqrt((x_target(irec)-dble(xstore(iglob)))**2 &
+(y_target(irec)-dble(ystore(iglob)))**2 &
+(z_target(irec)-dble(zstore(iglob)))**2)
@@ -291,16 +340,185 @@
ix_initial_guess(irec) = i
iy_initial_guess(irec) = j
iz_initial_guess(irec) = k
+
+ xi_receiver(irec) = dble(ix_initial_guess(irec))
+ eta_receiver(irec) = dble(iy_initial_guess(irec))
+ gamma_receiver(irec) = dble(iz_initial_guess(irec))
+ x_found(irec) = xstore(iglob)
+ y_found(irec) = ystore(iglob)
+ z_found(irec) = zstore(iglob)
endif
enddo
enddo
enddo
+
+! compute final distance between asked and found (converted to km)
+ final_distance(irec) = dsqrt((x_target(irec)-x_found(irec))**2 + &
+ (y_target(irec)-y_found(irec))**2 + (z_target(irec)-z_found(irec))**2)
! endif
! end of loop on all the spectral elements in current slice
enddo
+ if (ispec_selected_rec(irec) == 0) then
+ final_distance(irec) = HUGEVAL
+ endif
+
+! get normal to the face of the hexaedra if receiver is on the surface
+ if (USE_EXTERNAL_MESH .and. (.not. RECVS_CAN_BE_BURIED_EXT_MESH) .and. &
+ .not. (ispec_selected_rec(irec) == 0)) then
+ pt0_ix = -1
+ pt0_iy = -1
+ pt0_iz = -1
+ pt1_ix = -1
+ pt1_iy = -1
+ pt1_iz = -1
+ pt2_ix = -1
+ pt2_iy = -1
+ pt2_iz = -1
+! we get two vectors of the face (three points) to compute the normal
+ if (ix_initial_guess(irec) == 1 .and. &
+ iglob_is_surface_external_mesh(ibool(1,2,2,ispec_selected_rec(irec)))) then
+ pt0_ix = 1
+ pt0_iy = NGLLY
+ pt0_iz = 1
+ pt1_ix = 1
+ pt1_iy = 1
+ pt1_iz = 1
+ pt2_ix = 1
+ pt2_iy = NGLLY
+ pt2_iz = NGLLZ
+ endif
+ if (ix_initial_guess(irec) == NGLLX .and. &
+ iglob_is_surface_external_mesh(ibool(NGLLX,2,2,ispec_selected_rec(irec)))) then
+ pt0_ix = NGLLX
+ pt0_iy = 1
+ pt0_iz = 1
+ pt1_ix = NGLLX
+ pt1_iy = NGLLY
+ pt1_iz = 1
+ pt2_ix = NGLLX
+ pt2_iy = 1
+ pt2_iz = NGLLZ
+ endif
+ if (iy_initial_guess(irec) == 1 .and. &
+ iglob_is_surface_external_mesh(ibool(2,1,2,ispec_selected_rec(irec)))) then
+ pt0_ix = 1
+ pt0_iy = 1
+ pt0_iz = 1
+ pt1_ix = NGLLX
+ pt1_iy = 1
+ pt1_iz = 1
+ pt2_ix = 1
+ pt2_iy = 1
+ pt2_iz = NGLLZ
+ endif
+ if (iy_initial_guess(irec) == NGLLY .and. &
+ iglob_is_surface_external_mesh(ibool(2,NGLLY,2,ispec_selected_rec(irec)))) then
+ pt0_ix = NGLLX
+ pt0_iy = NGLLY
+ pt0_iz = 1
+ pt1_ix = 1
+ pt1_iy = NGLLY
+ pt1_iz = 1
+ pt2_ix = NGLLX
+ pt2_iy = NGLLY
+ pt2_iz = NGLLZ
+ endif
+ if (iz_initial_guess(irec) == 1 .and. &
+ iglob_is_surface_external_mesh(ibool(2,2,1,ispec_selected_rec(irec)))) then
+ pt0_ix = NGLLX
+ pt0_iy = 1
+ pt0_iz = 1
+ pt1_ix = 1
+ pt1_iy = 1
+ pt1_iz = 1
+ pt2_ix = NGLLX
+ pt2_iy = NGLLY
+ pt2_iz = 1
+ endif
+ if (iz_initial_guess(irec) == NGLLZ .and. &
+ iglob_is_surface_external_mesh(ibool(2,2,NGLLZ,ispec_selected_rec(irec)))) then
+ pt0_ix = 1
+ pt0_iy = 1
+ pt0_iz = NGLLZ
+ pt1_ix = NGLLX
+ pt1_iy = 1
+ pt1_iz = NGLLZ
+ pt2_ix = 1
+ pt2_iy = NGLLY
+ pt2_iz = NGLLZ
+ endif
+
+ if (pt0_ix<0 .or.pt0_iy<0 .or. pt0_iz<0 .or. &
+ pt1_ix<0 .or. pt1_iy<0 .or. pt1_iz<0 .or. &
+ pt2_ix<0 .or. pt2_iy<0 .or. pt2_iz<0) then
+ stop 'error in computing normal for receivers.'
+ endif
+
+ u_vector(1) = xstore(ibool(pt1_ix,pt1_iy,pt1_iz,ispec_selected_rec(irec))) &
+ - xstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+ u_vector(2) = ystore(ibool(pt1_ix,pt1_iy,pt1_iz,ispec_selected_rec(irec))) &
+ - ystore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+ u_vector(3) = zstore(ibool(pt1_ix,pt1_iy,pt1_iz,ispec_selected_rec(irec))) &
+ - zstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+ v_vector(1) = xstore(ibool(pt2_ix,pt2_iy,pt2_iz,ispec_selected_rec(irec))) &
+ - xstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+ v_vector(2) = ystore(ibool(pt2_ix,pt2_iy,pt2_iz,ispec_selected_rec(irec))) &
+ - ystore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+ v_vector(3) = zstore(ibool(pt2_ix,pt2_iy,pt2_iz,ispec_selected_rec(irec))) &
+ - zstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_rec(irec)))
+
+! cross product
+ w_vector(1) = u_vector(2)*v_vector(3) - u_vector(3)*v_vector(2)
+ w_vector(2) = u_vector(3)*v_vector(1) - u_vector(1)*v_vector(3)
+ w_vector(3) = u_vector(1)*v_vector(2) - u_vector(2)*v_vector(1)
+
+! normalize vector w
+ w_vector(:) = w_vector(:)/sqrt(w_vector(1)**2+w_vector(2)**2+w_vector(3)**2)
+
+! build the two other vectors for a direct base : we normalize u, and v=w^u
+ u_vector(:) = u_vector(:)/sqrt(u_vector(1)**2+u_vector(2)**2+u_vector(3)**2)
+ v_vector(1) = w_vector(2)*u_vector(3) - w_vector(3)*u_vector(2)
+ v_vector(2) = w_vector(3)*u_vector(1) - w_vector(1)*u_vector(3)
+ v_vector(3) = w_vector(1)*u_vector(2) - w_vector(2)*u_vector(1)
+
+! build rotation matrice nu for seismograms
+ if (EXT_MESH_RECV_NORMAL) then
+! East (u)
+ nu(1,1,irec) = u_vector(1)
+ nu(1,2,irec) = v_vector(1)
+ nu(1,3,irec) = w_vector(1)
+
+! North (v)
+ nu(2,1,irec) = u_vector(2)
+ nu(2,2,irec) = v_vector(2)
+ nu(2,3,irec) = w_vector(2)
+
+! Vertical (w)
+ nu(3,1,irec) = u_vector(3)
+ nu(3,2,irec) = v_vector(3)
+ nu(3,3,irec) = w_vector(3)
+ else
+! East
+ nu(1,1,irec) = 1.d0
+ nu(1,2,irec) = 0.d0
+ nu(1,3,irec) = 0.d0
+
+! North
+ nu(2,1,irec) = 0.d0
+ nu(2,2,irec) = 1.d0
+ nu(2,3,irec) = 0.d0
+
+! Vertical
+ nu(3,1,irec) = 0.d0
+ nu(3,2,irec) = 0.d0
+ nu(3,3,irec) = 1.d0
+ endif
+
+ endif ! of if (USE_EXTERNAL_MESH .and. (.not. RECVS_CAN_BE_BURIED_EXT_MESH))
+
! end of loop on all the stations
enddo
@@ -311,6 +529,8 @@
! find the best (xi,eta,gamma) for each receiver
! ****************************************
+ if(.not. FASTER_RECEIVERS_POINTS_ONLY) then
+
! loop on all the receivers to iterate in that slice
do irec = 1,nrec
@@ -424,6 +644,8 @@
enddo
+ endif ! of if (.not. FASTER_RECEIVERS_POINTS_ONLY)
+
! synchronize all the processes to make sure all the estimates are available
call sync_all()
@@ -437,6 +659,7 @@
call gather_all_dp(x_found,nrec,x_found_all,nrec,NPROC)
call gather_all_dp(y_found,nrec,y_found_all,nrec,NPROC)
call gather_all_dp(z_found,nrec,z_found_all,nrec,NPROC)
+ call gather_all_dp(nu,3*3*nrec,nu_all,3*3*nrec,NPROC)
! this is executed by main process only
if(myrank == 0) then
@@ -459,6 +682,7 @@
x_found(irec) = x_found_all(irec,iprocloop)
y_found(irec) = y_found_all(irec,iprocloop)
z_found(irec) = z_found_all(irec,iprocloop)
+ nu(:,:,irec) = nu_all(:,:,irec,iprocloop)
endif
enddo
final_distance(irec) = distmin
@@ -481,7 +705,15 @@
write(IMAIN,*) 'closest estimate found: ',sngl(final_distance(irec)),' m away'
write(IMAIN,*) ' in slice ',islice_selected_rec(irec),' in element ',ispec_selected_rec(irec)
- write(IMAIN,*) ' at xi,eta,gamma coordinates = ',xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec)
+ if(FASTER_RECEIVERS_POINTS_ONLY) then
+ write(IMAIN,*) 'in point i,j,k = ',nint(xi_receiver(irec)),nint(eta_receiver(irec)),nint(gamma_receiver(irec))
+ !write(IMAIN,*) 'in point i,j,k = ',x_found(irec),y_found(irec),z_found(irec)
+ write(IMAIN,*) 'nu1 = ',nu(1,:,irec)
+ write(IMAIN,*) 'nu2 = ',nu(2,:,irec)
+ write(IMAIN,*) 'nu3 = ',nu(3,:,irec)
+ else
+ write(IMAIN,*) ' at xi,eta,gamma coordinates = ',xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec)
+ endif
! add warning if estimate is poor
! (usually means receiver outside the mesh given by the user)
Modified: seismo/3D/SPECFEM3D/trunk/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/locate_source.f90 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/locate_source.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -36,7 +36,9 @@
LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX,Z_DEPTH_BLOCK, &
TOPOGRAPHY,itopo_bathy,UTM_PROJECTION_ZONE, &
PRINT_SOURCE_TIME_FUNCTION,SUPPRESS_UTM_PROJECTION, &
- NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO)
+ NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO, &
+ nu_source,iglob_is_surface_external_mesh,ispec_is_surface_external_mesh &
+ )
implicit none
@@ -68,7 +70,7 @@
integer iprocloop
- integer i,j,k,ispec,iglob,isource
+ integer i,j,k,ispec,iglob,isource,imin,imax,jmin,jmax,kmin,kmax
double precision, dimension(NSOURCES) :: utm_x_source,utm_y_source
double precision dist
@@ -112,11 +114,13 @@
integer, dimension(NGATHER_SOURCES,0:NPROC-1) :: ispec_selected_source_all
double precision, dimension(NGATHER_SOURCES,0:NPROC-1) :: xi_source_all,eta_source_all,gamma_source_all, &
final_distance_source_all,x_found_source_all,y_found_source_all,z_found_source_all
+ double precision, dimension(3,3,NGATHER_SOURCES,0:NPROC-1) :: nu_source_all
double precision hdur(NSOURCES), hdur_gaussian(NSOURCES), t0
double precision, dimension(NSOURCES) :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
+ double precision, dimension(3,3,NSOURCES) :: nu_source
integer icornerlong,icornerlat
double precision, dimension(NSOURCES) :: lat,long,depth,elevation
@@ -127,6 +131,12 @@
double precision, dimension(NSOURCES) :: x_found_source,y_found_source,z_found_source
double precision distmin
+! for surface locating and normal computing with external mesh
+ integer :: pt0_ix,pt0_iy,pt0_iz,pt1_ix,pt1_iy,pt1_iz,pt2_ix,pt2_iy,pt2_iz
+ real(kind=CUSTOM_REAL), dimension(3) :: u_vector,v_vector,w_vector
+ logical, dimension(NGLOB_AB) :: iglob_is_surface_external_mesh
+ logical, dimension(NSPEC_AB) :: ispec_is_surface_external_mesh
+
integer ix_initial_guess_source,iy_initial_guess_source,iz_initial_guess_source
! for calculation of source time function
@@ -157,13 +167,14 @@
! loop on all the sources
do isource = 1,NSOURCES
+ if (.not. USE_EXTERNAL_MESH) then
! check that the current source is inside the model
if(lat(isource) < LATITUDE_MIN .or. lat(isource) > LATITUDE_MAX .or. long(isource) < LONGITUDE_MIN &
.or. long(isource) > LONGITUDE_MAX) call exit_MPI(myrank,'the current source is outside the model')
if(depth(isource) >= dabs(Z_DEPTH_BLOCK/1000.d0)) &
call exit_MPI(myrank,'the current source is below the bottom of the model')
-
+ endif
!
! r -> z, theta -> -y, phi -> x
!
@@ -183,8 +194,27 @@
Mxy(isource) = - moment_tensor(6,isource)
call utm_geo(long(isource),lat(isource),utm_x_source(isource),utm_y_source(isource), &
- UTM_PROJECTION_ZONE,ILONGLAT2UTM,SUPPRESS_UTM_PROJECTION)
+ UTM_PROJECTION_ZONE,ILONGLAT2UTM,(SUPPRESS_UTM_PROJECTION .or. USE_EXTERNAL_MESH))
+! orientation consistent with the UTM projection
+
+! East
+ nu_source(1,1,isource) = 1.d0
+ nu_source(1,2,isource) = 0.d0
+ nu_source(1,3,isource) = 0.d0
+
+! North
+ nu_source(2,1,isource) = 0.d0
+ nu_source(2,2,isource) = 1.d0
+ nu_source(2,3,isource) = 0.d0
+
+! Vertical
+ nu_source(3,1,isource) = 0.d0
+ nu_source(3,2,isource) = 0.d0
+ nu_source(3,3,isource) = 1.d0
+
+ if (.not. USE_EXTERNAL_MESH) then
+
! compute elevation of topography at the epicenter
if(TOPOGRAPHY) then
@@ -230,20 +260,61 @@
z_target_source = - depth(isource)*1000.0d0 + elevation(isource)
if(myrank == 0) write(IOVTK,*) x_target_source,y_target_source,z_target_source
+ else
+
+ x_target_source = utm_x_source(isource)
+ y_target_source = utm_y_source(isource)
+ z_target_source = depth(isource)
+ if (myrank == 0) write(IOVTK,*) x_target_source, y_target_source, z_target_source
+
+ endif ! of if (.not. USE_EXTERNAL_MESH)
+
! set distance to huge initial value
distmin = HUGEVAL
+ ispec_selected_source(isource) = 0
+
do ispec=1,NSPEC_AB
+
+! define the interval in which we look for points
+ if(FASTER_SOURCES_POINTS_ONLY .and. USE_EXTERNAL_MESH) then
+ imin = 1
+ imax = NGLLX
+
+ jmin = 1
+ jmax = NGLLY
+
+ kmin = 1
+ kmax = NGLLZ
+
+ else
! loop only on points inside the element
! exclude edges to ensure this point is not shared with other elements
- do k=2,NGLLZ-1
- do j=2,NGLLY-1
- do i=2,NGLLX-1
+ imin = 2
+ imax = NGLLX - 1
-! keep this point if it is closer to the receiver
- iglob = ibool(i,j,k,ispec)
- dist=dsqrt((x_target_source-dble(xstore(iglob)))**2 &
+ jmin = 2
+ jmax = NGLLY - 1
+
+ kmin = 2
+ kmax = NGLLZ - 1
+ endif
+
+ do k = kmin,kmax
+ do j = jmin,jmax
+ do i = imin,imax
+
+ iglob = ibool(i,j,k,ispec)
+
+ if (USE_EXTERNAL_MESH .and. (.not. SOURCES_CAN_BE_BURIED_EXT_MESH)) then
+ if ((.not. iglob_is_surface_external_mesh(iglob)) .or. (.not. ispec_is_surface_external_mesh(ispec))) then
+ cycle
+ endif
+ endif
+
+! keep this point if it is closer to the source
+ dist=dsqrt((x_target_source-dble(xstore(iglob)))**2 &
+(y_target_source-dble(ystore(iglob)))**2 &
+(z_target_source-dble(zstore(iglob)))**2)
if(dist < distmin) then
@@ -252,6 +323,19 @@
ix_initial_guess_source = i
iy_initial_guess_source = j
iz_initial_guess_source = k
+
+! store xi,eta,gamma and x,y,z of point found
+ xi_source(isource) = dble(ix_initial_guess_source)
+ eta_source(isource) = dble(iy_initial_guess_source)
+ gamma_source(isource) = dble(iz_initial_guess_source)
+ x_found_source(isource) = xstore(iglob)
+ y_found_source(isource) = ystore(iglob)
+ z_found_source(isource) = zstore(iglob)
+
+! compute final distance between asked and found (converted to km)
+ final_distance_source(isource) = dsqrt((x_target_source-x_found_source(isource))**2 + &
+ (y_target_source-y_found_source(isource))**2 + (z_target_source-z_found_source(isource))**2)
+
endif
enddo
@@ -261,10 +345,153 @@
! end of loop on all the elements in current slice
enddo
+ if (ispec_selected_source(isource) == 0) then
+ final_distance_source(isource) = HUGEVAL
+ endif
+
+! get normal to the face of the hexaedra if receiver is on the surface
+ if (USE_EXTERNAL_MESH .and. (.not. SOURCES_CAN_BE_BURIED_EXT_MESH) .and. &
+ .not. (ispec_selected_source(isource) == 0)) then
+ pt0_ix = -1
+ pt0_iy = -1
+ pt0_iz = -1
+ pt1_ix = -1
+ pt1_iy = -1
+ pt1_iz = -1
+ pt2_ix = -1
+ pt2_iy = -1
+ pt2_iz = -1
+! we get two vectors of the face (three points) to compute the normal
+ if (xi_source(isource) == 1 .and. &
+ iglob_is_surface_external_mesh(ibool(1,2,2,ispec_selected_source(isource)))) then
+ pt0_ix = 1
+ pt0_iy = NGLLY
+ pt0_iz = 1
+ pt1_ix = 1
+ pt1_iy = 1
+ pt1_iz = 1
+ pt2_ix = 1
+ pt2_iy = NGLLY
+ pt2_iz = NGLLZ
+ endif
+ if (xi_source(isource) == NGLLX .and. &
+ iglob_is_surface_external_mesh(ibool(NGLLX,2,2,ispec_selected_source(isource)))) then
+ pt0_ix = NGLLX
+ pt0_iy = 1
+ pt0_iz = 1
+ pt1_ix = NGLLX
+ pt1_iy = NGLLY
+ pt1_iz = 1
+ pt2_ix = NGLLX
+ pt2_iy = 1
+ pt2_iz = NGLLZ
+ endif
+ if (eta_source(isource) == 1 .and. &
+ iglob_is_surface_external_mesh(ibool(2,1,2,ispec_selected_source(isource)))) then
+ pt0_ix = 1
+ pt0_iy = 1
+ pt0_iz = 1
+ pt1_ix = NGLLX
+ pt1_iy = 1
+ pt1_iz = 1
+ pt2_ix = 1
+ pt2_iy = 1
+ pt2_iz = NGLLZ
+ endif
+ if (eta_source(isource) == NGLLY .and. &
+ iglob_is_surface_external_mesh(ibool(2,NGLLY,2,ispec_selected_source(isource)))) then
+ pt0_ix = NGLLX
+ pt0_iy = NGLLY
+ pt0_iz = 1
+ pt1_ix = 1
+ pt1_iy = NGLLY
+ pt1_iz = 1
+ pt2_ix = NGLLX
+ pt2_iy = NGLLY
+ pt2_iz = NGLLZ
+ endif
+ if (gamma_source(isource) == 1 .and. &
+ iglob_is_surface_external_mesh(ibool(2,2,1,ispec_selected_source(isource)))) then
+ pt0_ix = NGLLX
+ pt0_iy = 1
+ pt0_iz = 1
+ pt1_ix = 1
+ pt1_iy = 1
+ pt1_iz = 1
+ pt2_ix = NGLLX
+ pt2_iy = NGLLY
+ pt2_iz = 1
+ endif
+ if (gamma_source(isource) == NGLLZ .and. &
+ iglob_is_surface_external_mesh(ibool(2,2,NGLLZ,ispec_selected_source(isource)))) then
+ pt0_ix = 1
+ pt0_iy = 1
+ pt0_iz = NGLLZ
+ pt1_ix = NGLLX
+ pt1_iy = 1
+ pt1_iz = NGLLZ
+ pt2_ix = 1
+ pt2_iy = NGLLY
+ pt2_iz = NGLLZ
+ endif
+
+ if (pt0_ix<0 .or.pt0_iy<0 .or. pt0_iz<0 .or. &
+ pt1_ix<0 .or. pt1_iy<0 .or. pt1_iz<0 .or. &
+ pt2_ix<0 .or. pt2_iy<0 .or. pt2_iz<0) then
+ stop 'error in computing normal for sources.'
+ endif
+
+ u_vector(1) = xstore(ibool(pt1_ix,pt1_iy,pt1_iz,ispec_selected_source(isource))) &
+ - xstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_source(isource)))
+ u_vector(2) = ystore(ibool(pt1_ix,pt1_iy,pt1_iz,ispec_selected_source(isource))) &
+ - ystore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_source(isource)))
+ u_vector(3) = zstore(ibool(pt1_ix,pt1_iy,pt1_iz,ispec_selected_source(isource))) &
+ - zstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_source(isource)))
+ v_vector(1) = xstore(ibool(pt2_ix,pt2_iy,pt2_iz,ispec_selected_source(isource))) &
+ - xstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_source(isource)))
+ v_vector(2) = ystore(ibool(pt2_ix,pt2_iy,pt2_iz,ispec_selected_source(isource))) &
+ - ystore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_source(isource)))
+ v_vector(3) = zstore(ibool(pt2_ix,pt2_iy,pt2_iz,ispec_selected_source(isource))) &
+ - zstore(ibool(pt0_ix,pt0_iy,pt0_iz,ispec_selected_source(isource)))
+
+! cross product
+ w_vector(1) = u_vector(2)*v_vector(3) - u_vector(3)*v_vector(2)
+ w_vector(2) = u_vector(3)*v_vector(1) - u_vector(1)*v_vector(3)
+ w_vector(3) = u_vector(1)*v_vector(2) - u_vector(2)*v_vector(1)
+
+! normalize vector w
+ w_vector(:) = w_vector(:)/sqrt(w_vector(1)**2+w_vector(2)**2+w_vector(3)**2)
+
+! build the two other vectors for a disourcet base : we normalize u, and v=w^u
+ u_vector(:) = u_vector(:)/sqrt(u_vector(1)**2+u_vector(2)**2+u_vector(3)**2)
+ v_vector(1) = w_vector(2)*u_vector(3) - w_vector(3)*u_vector(2)
+ v_vector(2) = w_vector(3)*u_vector(1) - w_vector(1)*u_vector(3)
+ v_vector(3) = w_vector(1)*u_vector(2) - w_vector(2)*u_vector(1)
+
+! build rotation matrice nu for seismograms
+! East (u)
+ nu_source(1,1,isource) = u_vector(1)
+ nu_source(1,2,isource) = v_vector(1)
+ nu_source(1,3,isource) = w_vector(1)
+
+! North (v)
+ nu_source(2,1,isource) = u_vector(2)
+ nu_source(2,2,isource) = v_vector(2)
+ nu_source(2,3,isource) = w_vector(2)
+
+! Vertical (w)
+ nu_source(3,1,isource) = u_vector(3)
+ nu_source(3,2,isource) = v_vector(3)
+ nu_source(3,3,isource) = w_vector(3)
+
+ endif ! of if (USE_EXTERNAL_MESH .and. (.not. RECEIVERS_CAN_BE_BURIED_EXT_MESH))
+
! *******************************************
! find the best (xi,eta,gamma) for the source
! *******************************************
+ if(.not. FASTER_SOURCES_POINTS_ONLY) then
+
! use initial guess in xi, eta and gamma
xi = xigll(ix_initial_guess_source)
eta = yigll(iy_initial_guess_source)
@@ -360,6 +587,8 @@
final_distance_source(isource) = dsqrt((x_target_source-x_found_source(isource))**2 + &
(y_target_source-y_found_source(isource))**2 + (z_target_source-z_found_source(isource))**2)
+ endif ! of if (.not. FASTER_SOURCES_POINTS_ONLY)
+
! end of loop on all the sources
enddo
@@ -373,7 +602,7 @@
ispec_selected_source_all(:,:) = -1
- call gather_all_i(ispec_selected_source(ns:ne),ng,ispec_selected_source_all(1:ng,:),ng)
+ call gather_all_i(ispec_selected_source(ns:ne),ng,ispec_selected_source_all(1:ng,:),ng,NPROC)
call gather_all_dp(xi_source(ns:ne),ng,xi_source_all(1:ng,:),ng,NPROC)
call gather_all_dp(eta_source(ns:ne),ng,eta_source_all(1:ng,:),ng,NPROC)
@@ -382,6 +611,7 @@
call gather_all_dp(x_found_source(ns:ne),ng,x_found_source_all(1:ng,:),ng,NPROC)
call gather_all_dp(y_found_source(ns:ne),ng,y_found_source_all(1:ng,:),ng,NPROC)
call gather_all_dp(z_found_source(ns:ne),ng,z_found_source_all(1:ng,:),ng,NPROC)
+ call gather_all_dp(nu_source(:,:,ns:ne),3*3*ng,nu_source_all(:,:,1:ng,:),3*3*ng,NPROC)
! this is executed by main process only
if(myrank == 0) then
@@ -406,6 +636,7 @@
x_found_source(isource) = x_found_source_all(is,iprocloop)
y_found_source(isource) = y_found_source_all(is,iprocloop)
z_found_source(isource) = z_found_source_all(is,iprocloop)
+ nu_source(:,:,isource) = nu_source_all(:,:,isource,iprocloop)
endif
enddo
final_distance_source(isource) = distmin
@@ -428,9 +659,19 @@
write(IMAIN,*) 'source located in slice ',islice_selected_source(isource)
write(IMAIN,*) ' in element ',ispec_selected_source(isource)
write(IMAIN,*)
- write(IMAIN,*) ' xi coordinate of source in that element: ',xi_source(isource)
- write(IMAIN,*) ' eta coordinate of source in that element: ',eta_source(isource)
- write(IMAIN,*) 'gamma coordinate of source in that element: ',gamma_source(isource)
+ if(FASTER_SOURCES_POINTS_ONLY) then
+ write(IMAIN,*) ' xi coordinate of source in that element: ',nint(xi_source(isource))
+ write(IMAIN,*) ' eta coordinate of source in that element: ',nint(eta_source(isource))
+ write(IMAIN,*) 'gamma coordinate of source in that element: ',nint(gamma_source(isource))
+ write(IMAIN,*) 'nu1 = ',nu_source(1,:,isource)
+ write(IMAIN,*) 'nu2 = ',nu_source(2,:,isource)
+ write(IMAIN,*) 'nu3 = ',nu_source(3,:,isource)
+ write(IMAIN,*) 'at (x,y,z) coordinates = ',x_found_source(isource),y_found_source(isource),z_found_source(isource)
+ else
+ write(IMAIN,*) ' xi coordinate of source in that element: ',xi_source(isource)
+ write(IMAIN,*) ' eta coordinate of source in that element: ',eta_source(isource)
+ write(IMAIN,*) 'gamma coordinate of source in that element: ',gamma_source(isource)
+ endif
! add message if source is a Heaviside
if(hdur(isource) < 5.*DT) then
Modified: seismo/3D/SPECFEM3D/trunk/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/meshfem3D.f90 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/meshfem3D.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -228,6 +228,25 @@
character(len=MAX_LENGTH_NETWORK_NAME) network_name
character(len=150) rec_filename,filtered_rec_filename,dummystring
+! for Databases of external meshes
+ character(len=150) prname
+ integer :: dummy_node
+ integer :: dummy_elmnt
+ integer :: ispec, inode, num_interface, ie
+ integer :: nnodes_ext_mesh, nelmnts_ext_mesh
+ integer :: ninterface_ext_mesh
+ integer :: max_interface_size_ext_mesh
+ integer, dimension(:), allocatable :: my_neighbours_ext_mesh
+ integer, dimension(:), allocatable :: my_nelmnts_neighbours_ext_mesh
+ integer, dimension(:,:,:), allocatable :: my_interfaces_ext_mesh
+ integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh
+ integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh_dummy
+ integer, dimension(:), allocatable :: ibool_interface_ext_mesh_dummy
+ integer, dimension(:), allocatable :: nibool_interfaces_ext_mesh
+ double precision, dimension(:,:), allocatable :: nodes_coords_ext_mesh
+ integer, dimension(:,:), allocatable :: elmnts_ext_mesh
+ integer, dimension(:), allocatable :: mat_ext_mesh
+
! ************** PROGRAM STARTS HERE **************
! sizeprocs returns number of processes started (should be equal to NPROC).
@@ -282,9 +301,16 @@
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,USE_REGULAR_MESH)
+! info about external mesh simulation
+! nlegoff -- should be put in compute_parameters and read_parameter_file for clarity
+ if (USE_EXTERNAL_MESH) then
+ NPROC = sizeprocs
+ endif
+
! check that the code is running with the requested nb of processes
if(sizeprocs /= NPROC) call exit_MPI(myrank,'wrong number of MPI processes')
+ if (.not. USE_EXTERNAL_MESH) then
! dynamic allocation of mesh arrays
allocate(rns(0:2*NER))
@@ -330,6 +356,8 @@
write(IMAIN,'(a1)',advance='yes') ' '
enddo
endif
+
+ endif ! end of (.not. USE_EXTERNAL_MESH)
if(myrank == 0) then
write(IMAIN,*) 'This is process ',myrank
@@ -362,6 +390,7 @@
! for the number of standard linear solids for attenuation
if(N_SLS /= 3) call exit_MPI(myrank,'number of SLS must be 3')
+ if (.not. USE_EXTERNAL_MESH) then
! check that Poisson's ratio in Gocad block is fine
if(VP_VS_RATIO_GOCAD_TOP < sqrt(2.) .or. VP_VS_RATIO_GOCAD_BOTTOM < sqrt(2.))&
call exit_MPI(myrank,'vp/vs ratio in Gocad block is too small')
@@ -383,6 +412,8 @@
if(mod(NEX_ETA/8,NPROC_ETA) /= 0) call exit_MPI(myrank,'NEX_ETA must be a multiple of 8*NPROC_ETA for a non-regular mesh')
endif
+ endif ! end of (.not. USE_EXTERNAL_MESH)
+
if(myrank == 0) then
write(IMAIN,*) 'region selected:'
@@ -498,6 +529,8 @@
close(55)
endif
+ if (.not. USE_EXTERNAL_MESH) then
+
! get addressing for this process
iproc_xi = iproc_xi_slice(myrank)
iproc_eta = iproc_eta_slice(myrank)
@@ -672,6 +705,8 @@
enddo
enddo
+ endif ! end of (.not. USE_EXTERNAL_MESH)
+
if(myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) '**************************'
@@ -685,6 +720,46 @@
area_local_bottom = ZERO
area_local_top = ZERO
+! read databases about external mesh simulation
+! nlegoff --
+ if (USE_EXTERNAL_MESH) then
+ call create_name_database(prname,myrank,LOCAL_PATH)
+ open(unit=IIN,file=prname(1:len_trim(prname))//'Database',status='old',action='read',form='formatted')
+ read(IIN,*) nnodes_ext_mesh
+ allocate(nodes_coords_ext_mesh(NDIM,nnodes_ext_mesh))
+ do inode = 1, nnodes_ext_mesh
+ read(IIN,*) dummy_node, nodes_coords_ext_mesh(1,inode), nodes_coords_ext_mesh(2,inode), nodes_coords_ext_mesh(3,inode)
+ enddo
+
+ read(IIN,*) nelmnts_ext_mesh
+ allocate(elmnts_ext_mesh(esize,nelmnts_ext_mesh))
+ allocate(mat_ext_mesh(nelmnts_ext_mesh))
+ do ispec = 1, nelmnts_ext_mesh
+ read(IIN,*) dummy_elmnt, mat_ext_mesh(ispec), &
+ elmnts_ext_mesh(1,ispec), elmnts_ext_mesh(2,ispec), elmnts_ext_mesh(3,ispec), elmnts_ext_mesh(4,ispec), &
+ elmnts_ext_mesh(5,ispec), elmnts_ext_mesh(6,ispec), elmnts_ext_mesh(7,ispec), elmnts_ext_mesh(8,ispec)
+ enddo
+ NSPEC_AB = nelmnts_ext_mesh
+
+ read(IIN,*) ninterface_ext_mesh, max_interface_size_ext_mesh
+ allocate(my_neighbours_ext_mesh(ninterface_ext_mesh))
+ allocate(my_nelmnts_neighbours_ext_mesh(ninterface_ext_mesh))
+ allocate(my_interfaces_ext_mesh(6,max_interface_size_ext_mesh,ninterface_ext_mesh))
+ allocate(ibool_interfaces_ext_mesh(NGLLX*NGLLX*max_interface_size_ext_mesh,ninterface_ext_mesh))
+ allocate(nibool_interfaces_ext_mesh(ninterface_ext_mesh))
+ do num_interface = 1, ninterface_ext_mesh
+ read(IIN,*) my_neighbours_ext_mesh(num_interface), my_nelmnts_neighbours_ext_mesh(num_interface)
+ do ie = 1, my_nelmnts_neighbours_ext_mesh(num_interface)
+ read(IIN,*) my_interfaces_ext_mesh(1,ie,num_interface), my_interfaces_ext_mesh(2,ie,num_interface), &
+ my_interfaces_ext_mesh(3,ie,num_interface), my_interfaces_ext_mesh(4,ie,num_interface), &
+ my_interfaces_ext_mesh(5,ie,num_interface), my_interfaces_ext_mesh(6,ie,num_interface)
+ enddo
+ enddo
+
+ close(IIN)
+
+ endif
+
! assign theoretical number of elements
nspec = NSPEC_AB
@@ -705,6 +780,20 @@
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! create all the regions of the mesh
+ if (USE_EXTERNAL_MESH) then
+ call create_regions_mesh_ext_mesh(ibool, &
+ xstore,ystore,zstore,npx,npy,iproc_xi,iproc_eta,nspec, &
+ volume_local,area_local_bottom,area_local_top, &
+ NGLOB_AB,npointot, &
+ myrank,LOCAL_PATH, &
+ nnodes_ext_mesh,nelmnts_ext_mesh, &
+ nodes_coords_ext_mesh,elmnts_ext_mesh,mat_ext_mesh, &
+ ninterface_ext_mesh,max_interface_size_ext_mesh, &
+ my_neighbours_ext_mesh,my_nelmnts_neighbours_ext_mesh,my_interfaces_ext_mesh, &
+ ibool_interfaces_ext_mesh,nibool_interfaces_ext_mesh &
+ )
+
+ else
call create_regions_mesh(xgrid,ygrid,zgrid,ibool,idoubling, &
xstore,ystore,zstore,npx,npy, &
iproc_xi,iproc_eta,nspec, &
@@ -722,7 +811,7 @@
IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,MOHO_MAP_LUPEI, &
ANISOTROPY,SAVE_MESH_FILES,SUPPRESS_UTM_PROJECTION, &
ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO,NX_TOPO,NY_TOPO,USE_REGULAR_MESH)
-
+ endif
! print min and max of topography included
if(TOPOGRAPHY) then
@@ -831,7 +920,8 @@
open(unit=IIN,file=rec_filename,status='old',action='read')
do irec = 1,nrec
read(IIN,*) station_name,network_name,stlat,stlon,stele,stbur
- if(stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+ if((stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+ .or. USE_EXTERNAL_MESH) &
nrec_filtered = nrec_filtered + 1
enddo
close(IIN)
@@ -849,7 +939,8 @@
do irec = 1,nrec
read(IIN,*) station_name,network_name,stlat,stlon,stele,stbur
- if(stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+ if((stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+ .or. USE_EXTERNAL_MESH) &
write(IOUT,*) station_name(1:len_trim(station_name)),' ',network_name(1:len_trim(network_name)),' ', &
sngl(stlat),' ',sngl(stlon), ' ', sngl(stele), ' ', sngl(stbur)
enddo
Modified: seismo/3D/SPECFEM3D/trunk/parallel.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/parallel.f90 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/parallel.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -191,6 +191,33 @@
!----
!
+ subroutine gatherv_all_cr(sendbuf, sendcnt, recvbuf, recvcount, recvoffset,recvcounttot, NPROC)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+ include "precision.h"
+
+ integer sendcnt,recvcounttot,NPROC
+ integer, dimension(NPROC) :: recvcount,recvoffset
+ real(kind=CUSTOM_REAL), dimension(sendcnt) :: sendbuf
+ real(kind=CUSTOM_REAL), dimension(recvcounttot) :: recvbuf
+
+ integer ier
+
+ call MPI_GATHERV(sendbuf,sendcnt,CUSTOM_MPI_TYPE, &
+ recvbuf,recvcount,recvoffset,CUSTOM_MPI_TYPE, &
+ 0,MPI_COMM_WORLD,ier)
+
+ end subroutine gatherv_all_cr
+
+!
+!----
+!
+
subroutine init()
implicit none
@@ -401,3 +428,132 @@
proc_null = MPI_PROC_NULL
end function proc_null
+
+!
+!----
+!
+
+ subroutine issend_cr(sendbuf, sendcount, dest, sendtag, req)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+ include "precision.h"
+
+ integer sendcount, dest, sendtag, req
+ real(kind=CUSTOM_REAL), dimension(sendcount) :: sendbuf
+
+! MPI status of messages to be received
+ integer msg_status(MPI_STATUS_SIZE)
+
+ integer ier
+
+ call MPI_ISSEND(sendbuf(1),sendcount,CUSTOM_MPI_TYPE,dest,sendtag, &
+ MPI_COMM_WORLD,req,ier)
+
+ end subroutine issend_cr
+
+!
+!----
+!
+
+ subroutine irecv_cr(recvbuf, recvcount, dest, recvtag, req)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+ include "precision.h"
+
+ integer recvcount, dest, recvtag, req
+ real(kind=CUSTOM_REAL), dimension(recvcount) :: recvbuf
+
+! MPI status of messages to be received
+ integer msg_status(MPI_STATUS_SIZE)
+
+ integer ier
+
+ call MPI_IRECV(recvbuf(1),recvcount,CUSTOM_MPI_TYPE,dest,recvtag, &
+ MPI_COMM_WORLD,req,ier)
+
+ end subroutine irecv_cr
+
+!
+!----
+!
+
+ subroutine issend_i(sendbuf, sendcount, dest, sendtag, req)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+ include "precision.h"
+
+ integer sendcount, dest, sendtag, req
+ integer, dimension(sendcount) :: sendbuf
+
+! MPI status of messages to be received
+ integer msg_status(MPI_STATUS_SIZE)
+
+ integer ier
+
+ call MPI_ISSEND(sendbuf(1),sendcount,MPI_INTEGER,dest,sendtag, &
+ MPI_COMM_WORLD,req,ier)
+
+ end subroutine issend_i
+
+!
+!----
+!
+
+ subroutine irecv_i(recvbuf, recvcount, dest, recvtag, req)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+ include "precision.h"
+
+ integer recvcount, dest, recvtag, req
+ integer, dimension(recvcount) :: recvbuf
+
+! MPI status of messages to be received
+ integer msg_status(MPI_STATUS_SIZE)
+
+ integer ier
+
+ call MPI_IRECV(recvbuf(1),recvcount,MPI_INTEGER,dest,recvtag, &
+ MPI_COMM_WORLD,req,ier)
+
+ end subroutine irecv_i
+
+!
+!----
+!
+
+ subroutine wait_req(req)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ integer :: req
+
+ integer, dimension(MPI_STATUS_SIZE) :: req_mpi_status
+
+ integer :: ier
+
+ call mpi_wait(req,req_mpi_status,ier)
+
+ end subroutine wait_req
Modified: seismo/3D/SPECFEM3D/trunk/read_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/read_arrays_solver.f90 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/read_arrays_solver.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -25,7 +25,7 @@
! read arrays created by the mesher
- subroutine read_arrays_solver(myrank,xstore,ystore,zstore, &
+ subroutine read_arrays_solver(myrank,NSPEC_AB,NGLOB_AB,xstore,ystore,zstore, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian, &
flag_sediments,not_fully_in_bedrock,rho_vp,rho_vs,ANISOTROPY, &
c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
@@ -41,41 +41,44 @@
integer myrank
+ integer NSPEC_AB
+ integer NGLOB_AB
+
logical OCEANS
character(len=150) LOCAL_PATH
! coordinates in single precision
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB_VAL) :: xstore,ystore,zstore
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian
logical ANISOTROPY
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO) :: &
c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store
! material properties
- real(kind=CUSTOM_REAL) kappastore(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL)
- real(kind=CUSTOM_REAL) mustore(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL)
+ real(kind=CUSTOM_REAL) kappastore(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
+ real(kind=CUSTOM_REAL) mustore(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
! flag for sediments
- logical not_fully_in_bedrock(NSPEC_AB_VAL)
- logical flag_sediments(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL)
+ logical not_fully_in_bedrock(NSPEC_AB)
+ logical flag_sediments(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
! Stacey
- real(kind=CUSTOM_REAL) rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL)
- real(kind=CUSTOM_REAL) rho_vs(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL)
+ real(kind=CUSTOM_REAL) rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
+ real(kind=CUSTOM_REAL) rho_vs(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
! mass matrix and additional ocean load mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB_VAL) :: rmass,rmass_ocean_load
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: rmass,rmass_ocean_load
! global addressing
- integer ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL)
+ integer ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB)
- integer idoubling(NSPEC_AB_VAL)
+ integer idoubling(NSPEC_AB)
! processor identification
character(len=150) prname
Copied: seismo/3D/SPECFEM3D/trunk/sort_array_coordinates.f90 (from rev 13200, seismo/3D/SPECFEM3D/branches/update_temporary/sort_array_coordinates.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/sort_array_coordinates.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/sort_array_coordinates.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -0,0 +1,235 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! subroutines to sort MPI buffers to assemble between chunks
+
+ subroutine sort_array_coordinates(npointot,x,y,z,ibool,iglob,loc,ifseg,nglob,ind,ninseg,iwork,work)
+
+! this routine MUST be in double precision to avoid sensitivity
+! to roundoff errors in the coordinates of the points
+
+ implicit none
+
+ include "constants.h"
+
+ integer npointot,nglob
+
+ integer ibool(npointot),iglob(npointot),loc(npointot)
+ integer ind(npointot),ninseg(npointot)
+ logical ifseg(npointot)
+ double precision x(npointot),y(npointot),z(npointot)
+ integer iwork(npointot)
+ double precision work(npointot)
+
+ integer ipoin,i,j
+ integer nseg,ioff,iseg,ig
+ double precision xtol
+
+! establish initial pointers
+ do ipoin=1,npointot
+ loc(ipoin)=ipoin
+ enddo
+
+! define a tolerance, normalized radius is 1., so let's use a small value
+ xtol = SMALLVAL_TOL
+
+ ifseg(:)=.false.
+
+ nseg=1
+ ifseg(1)=.true.
+ ninseg(1)=npointot
+
+ do j=1,NDIM
+
+! sort within each segment
+ ioff=1
+ do iseg=1,nseg
+ if(j == 1) then
+
+ call rank_buffers(x(ioff),ind,ninseg(iseg))
+
+ else if(j == 2) then
+
+ call rank_buffers(y(ioff),ind,ninseg(iseg))
+
+ else
+
+ call rank_buffers(z(ioff),ind,ninseg(iseg))
+
+ endif
+
+ call swap_all_buffers(ibool(ioff),loc(ioff), &
+ x(ioff),y(ioff),z(ioff),iwork,work,ind,ninseg(iseg))
+
+ ioff=ioff+ninseg(iseg)
+ enddo
+
+! check for jumps in current coordinate
+ if(j == 1) then
+ do i=2,npointot
+ if(dabs(x(i)-x(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ else if(j == 2) then
+ do i=2,npointot
+ if(dabs(y(i)-y(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ else
+ do i=2,npointot
+ if(dabs(z(i)-z(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ endif
+
+! count up number of different segments
+ nseg=0
+ do i=1,npointot
+ if(ifseg(i)) then
+ nseg=nseg+1
+ ninseg(nseg)=1
+ else
+ ninseg(nseg)=ninseg(nseg)+1
+ endif
+ enddo
+ enddo
+
+! assign global node numbers (now sorted lexicographically)
+ ig=0
+ do i=1,npointot
+ if(ifseg(i)) ig=ig+1
+ iglob(loc(i))=ig
+ enddo
+
+ nglob=ig
+
+ end subroutine sort_array_coordinates
+
+! -------------------- library for sorting routine ------------------
+
+! sorting routines put here in same file to allow for inlining
+
+ subroutine rank_buffers(A,IND,N)
+!
+! Use Heap Sort (Numerical Recipes)
+!
+ implicit none
+
+ integer n
+ double precision A(n)
+ integer IND(n)
+
+ integer i,j,l,ir,indx
+ double precision q
+
+ do j=1,n
+ IND(j)=j
+ enddo
+
+ if(n == 1) return
+
+ L=n/2+1
+ ir=n
+ 100 CONTINUE
+ IF(l>1) THEN
+ l=l-1
+ indx=ind(l)
+ q=a(indx)
+ ELSE
+ indx=ind(ir)
+ q=a(indx)
+ ind(ir)=ind(1)
+ ir=ir-1
+ if (ir == 1) then
+ ind(1)=indx
+ return
+ endif
+ ENDIF
+ i=l
+ j=l+l
+ 200 CONTINUE
+ IF(J <= IR) THEN
+ IF(J < IR) THEN
+ IF(A(IND(j)) < A(IND(j+1))) j=j+1
+ ENDIF
+ IF (q < A(IND(j))) THEN
+ IND(I)=IND(J)
+ I=J
+ J=J+J
+ ELSE
+ J=IR+1
+ ENDIF
+ goto 200
+ ENDIF
+ IND(I)=INDX
+ goto 100
+ end subroutine rank_buffers
+
+! -------------------------------------------------------------------
+
+ subroutine swap_all_buffers(IA,IB,A,B,C,IW,W,ind,n)
+!
+! swap arrays IA, IB, A, B and C according to addressing in array IND
+!
+ implicit none
+
+ integer n
+
+ integer IND(n)
+ integer IA(n),IB(n),IW(n)
+ double precision A(n),B(n),C(n),W(n)
+
+ integer i
+
+ do i=1,n
+ W(i)=A(i)
+ IW(i)=IA(i)
+ enddo
+
+ do i=1,n
+ A(i)=W(ind(i))
+ IA(i)=IW(ind(i))
+ enddo
+
+ do i=1,n
+ W(i)=B(i)
+ IW(i)=IB(i)
+ enddo
+
+ do i=1,n
+ B(i)=W(ind(i))
+ IB(i)=IW(ind(i))
+ enddo
+
+ do i=1,n
+ W(i)=C(i)
+ enddo
+
+ do i=1,n
+ C(i)=W(ind(i))
+ enddo
+
+ end subroutine swap_all_buffers
+
+
Modified: seismo/3D/SPECFEM3D/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/specfem3D.f90 2008-10-31 18:17:24 UTC (rev 13205)
+++ seismo/3D/SPECFEM3D/trunk/specfem3D.f90 2008-10-31 19:04:49 UTC (rev 13206)
@@ -185,13 +185,13 @@
! -----------------
! mesh parameters
- integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL) :: ibool
+ integer, dimension(:,:,:,:), allocatable :: ibool
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL) :: &
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB_VAL) :: xstore,ystore,zstore
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: xstore,ystore,zstore
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL) :: &
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
kappastore,mustore
! material properties in case of a fully anisotropic material
@@ -201,26 +201,26 @@
c36store,c44store,c45store,c46store,c55store,c56store,c66store
! flag for sediments
- logical not_fully_in_bedrock(NSPEC_AB_VAL)
- logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL) :: flag_sediments
+ logical, dimension(:), allocatable :: not_fully_in_bedrock
+ logical, dimension(:,:,:,:), allocatable :: flag_sediments
! Stacey
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL) :: rho_vp,rho_vs
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_vp,rho_vs
! local to global mapping
- integer, dimension(NSPEC_AB_VAL) :: idoubling
+ integer, dimension(:), allocatable :: idoubling
! mass matrix
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB_VAL) :: rmass
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass
! additional mass matrix for ocean load
! ocean load mass matrix is always allocated statically even if no oceans
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB_VAL) :: rmass_ocean_load
- logical, dimension(NGLOB_AB_VAL) :: updated_dof_ocean_load
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_ocean_load
+ logical, dimension(:), allocatable :: updated_dof_ocean_load
real(kind=CUSTOM_REAL) additional_term,force_normal_comp
! displacement, velocity, acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB_VAL) :: displ,veloc,accel
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: displ,veloc,accel
real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
@@ -301,6 +301,7 @@
real(kind=CUSTOM_REAL) stf_used
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: sourcearray
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: sourcearrays
+ double precision, dimension(:,:,:), allocatable :: nu_source
!ADJOINT
character(len=150) adj_source_file
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: adj_sourcearray
@@ -312,7 +313,7 @@
double precision, dimension(:), allocatable :: t_cmt,hdur,hdur_gaussian
double precision, dimension(:), allocatable :: utm_x_source,utm_y_source
double precision, external :: comp_source_time_function
- double precision :: t0
+ double precision :: t0,f0
! receiver information
character(len=150) rec_filename,filtered_rec_filename,dummystring
@@ -439,6 +440,51 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dvxdxl,dvxdyl,dvxdzl,dvydxl,dvydyl,dvydzl,dvzdxl,dvzdyl,dvzdzl
real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable:: div, curl_x, curl_y, curl_z
+! for assembling in case of external mesh
+ integer :: ninterfaces_ext_mesh
+ integer :: max_nibool_interfaces_ext_mesh
+ integer, dimension(:), allocatable :: my_neighbours_ext_mesh
+ integer, dimension(:), allocatable :: nibool_interfaces_ext_mesh
+ integer, dimension(:,:), allocatable :: ibool_interfaces_ext_mesh
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: buffer_send_vector_ext_mesh
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: buffer_recv_vector_ext_mesh
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: buffer_send_scalar_ext_mesh
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: buffer_recv_scalar_ext_mesh
+ integer, dimension(:), allocatable :: request_send_scalar_ext_mesh
+ integer, dimension(:), allocatable :: request_recv_scalar_ext_mesh
+ integer, dimension(:), allocatable :: request_send_vector_ext_mesh
+ integer, dimension(:), allocatable :: request_recv_vector_ext_mesh
+
+! for detecting surface receivers and source in case of external mesh
+ integer, dimension(:), allocatable :: valence_external_mesh
+ logical, dimension(:), allocatable :: iglob_is_surface_external_mesh
+ logical, dimension(:), allocatable :: ispec_is_surface_external_mesh
+ integer, dimension(:,:), allocatable :: buffer_send_scalar_i_ext_mesh
+ integer, dimension(:,:), allocatable :: buffer_recv_scalar_i_ext_mesh
+ integer :: nfaces_surface_external_mesh
+ integer :: nfaces_surface_glob_ext_mesh
+ integer,dimension(:),allocatable :: nfaces_perproc_surface_ext_mesh
+ integer,dimension(:),allocatable :: faces_surface_offset_ext_mesh
+ integer,dimension(:,:),allocatable :: faces_surface_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_x_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_y_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_z_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_x_all_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_y_all_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_z_all_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_ux_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_uy_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_uz_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_ux_all_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_uy_all_external_mesh
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_uz_all_external_mesh
+ integer :: ii,jj,kk
+
+! for communications overlapping
+ logical, dimension(:), allocatable :: ispec_is_inner_ext_mesh
+ logical, dimension(:), allocatable :: iglob_is_inner_ext_mesh
+ integer :: iinterface
+
! ************** PROGRAM STARTS HERE **************
! sizeprocs returns number of processes started
@@ -487,6 +533,19 @@
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
+! info about external mesh simulation
+! nlegoff -- should be put in compute_parameters and read_parameter_file for clarity
+ if (USE_EXTERNAL_MESH) then
+ NPROC = sizeprocs
+ DT = DT_ext_mesh
+ NSTEP = NSTEP_ext_mesh
+ call create_name_database(prname,myrank,LOCAL_PATH)
+ open(unit=27,file=prname(1:len_trim(prname))//'external_mesh.bin',status='old',action='read',form='unformatted')
+ read(27) NSPEC_AB
+ read(27) NGLOB_AB
+ close(27)
+ endif
+
! open main output file, only written to by process 0
if(myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) &
open(unit=IMAIN,file=trim(OUTPUT_FILES)//'/output_solver.txt',status='unknown')
@@ -540,6 +599,9 @@
! check that we have at least one source
if(NSOURCES < 1) call exit_MPI(myrank,'need at least one source')
+! info on the addressing scheme that is not used in case of an external mesh simulation
+ if (.not. USE_EXTERNAL_MESH) then
+
! open file with global slice number addressing
if(myrank == 0) then
open(unit=IIN,file=trim(OUTPUT_FILES)//'/addressing.txt',status='old',action='read')
@@ -565,10 +627,109 @@
! define maximum size for message buffers
NPOIN2DMAX_XY = max(NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX)
+ endif ! end of (.not. USE_EXTERNAL_MESH)
+
! start reading the databases
+! allocate arrays for storing the databases
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(xix(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(xiy(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(xiz(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(etax(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(etay(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(etaz(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(gammax(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(gammay(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(gammaz(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(xstore(NGLOB_AB))
+ allocate(ystore(NGLOB_AB))
+ allocate(zstore(NGLOB_AB))
+ allocate(kappastore(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(mustore(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(not_fully_in_bedrock(NSPEC_AB))
+ allocate(flag_sediments(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(rho_vs(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(idoubling(NSPEC_AB))
+ allocate(rmass(NGLOB_AB))
+ allocate(rmass_ocean_load(NGLOB_AB))
+ allocate(updated_dof_ocean_load(NGLOB_AB))
+ allocate(displ(NDIM,NGLOB_AB))
+ allocate(veloc(NDIM,NGLOB_AB))
+ allocate(accel(NDIM,NGLOB_AB))
+
+! info about external mesh simulation
+! nlegoff -- should be put in read_arrays_solver and read_arrays_buffer_solver for clarity
+ if (USE_EXTERNAL_MESH) then
+ call create_name_database(prname,myrank,LOCAL_PATH)
+ open(unit=27,file=prname(1:len_trim(prname))//'external_mesh.bin',status='old',action='read',form='unformatted')
+ read(27) NSPEC_AB
+ read(27) NGLOB_AB
+ read(27) xix
+ read(27) xiy
+ read(27) xiz
+ read(27) etax
+ read(27) etay
+ read(27) etaz
+ read(27) gammax
+ read(27) gammay
+ read(27) gammaz
+ read(27) jacobian
+ read(27) kappastore
+ read(27) mustore
+ read(27) rmass
+ read(27) ibool
+ read(27) xstore
+ read(27) ystore
+ read(27) zstore
+
+ read(27) ninterfaces_ext_mesh
+ read(27) max_nibool_interfaces_ext_mesh
+ allocate(my_neighbours_ext_mesh(ninterfaces_ext_mesh))
+ allocate(nibool_interfaces_ext_mesh(ninterfaces_ext_mesh))
+ allocate(ibool_interfaces_ext_mesh(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh))
+ read(27) my_neighbours_ext_mesh
+ read(27) nibool_interfaces_ext_mesh
+ read(27) ibool_interfaces_ext_mesh
+
+ allocate(buffer_send_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh))
+ allocate(buffer_recv_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh))
+ allocate(buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh))
+ allocate(buffer_recv_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh))
+ allocate(request_send_vector_ext_mesh(ninterfaces_ext_mesh))
+ allocate(request_recv_vector_ext_mesh(ninterfaces_ext_mesh))
+ allocate(request_send_scalar_ext_mesh(ninterfaces_ext_mesh))
+ allocate(request_recv_scalar_ext_mesh(ninterfaces_ext_mesh))
+ close(27)
+
+! locate inner and outer elements
+ allocate(ispec_is_inner_ext_mesh(NSPEC_AB))
+ allocate(iglob_is_inner_ext_mesh(NGLOB_AB))
+ ispec_is_inner_ext_mesh(:) = .true.
+ iglob_is_inner_ext_mesh(:) = .true.
+ do iinterface = 1, ninterfaces_ext_mesh
+ do i = 1, nibool_interfaces_ext_mesh(iinterface)
+ iglob = ibool_interfaces_ext_mesh(i,iinterface)
+ iglob_is_inner_ext_mesh(iglob) = .false.
+ enddo
+ enddo
+ do ispec = 1, NSPEC_AB
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+ ispec_is_inner_ext_mesh(ispec) = iglob_is_inner_ext_mesh(iglob) .and. ispec_is_inner_ext_mesh(ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+
+ else
+
! read arrays created by the mesher
- call read_arrays_solver(myrank,xstore,ystore,zstore, &
+ call read_arrays_solver(myrank,NSPEC_AB,NGLOB_AB,xstore,ystore,zstore, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian, &
flag_sediments,not_fully_in_bedrock,rho_vp,rho_vs,ANISOTROPY, &
c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
@@ -577,7 +738,7 @@
kappastore,mustore,ibool,idoubling,rmass,rmass_ocean_load,LOCAL_PATH,OCEANS)
! check that the number of points in this slice is correct
- if(minval(ibool(:,:,:,:)) /= 1 .or. maxval(ibool(:,:,:,:)) /= NGLOB_AB_VAL) &
+ if(minval(ibool(:,:,:,:)) /= 1 .or. maxval(ibool(:,:,:,:)) /= NGLOB_AB) &
call exit_MPI(myrank,'incorrect global numbering: iboolmax does not equal NGLOB')
! read 2-D addressing for summation between slices with MPI
@@ -586,6 +747,8 @@
npoin2D_xi,npoin2D_eta, &
NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,LOCAL_PATH)
+ endif ! end of (USE_EXTERNAL_MESH)
+
! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
if(myrank == 0) then
@@ -619,6 +782,7 @@
if (ATTENUATION .and. ((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3)) &
call create_name_database(prname_Q,myrank,LOCAL_PATH_Q)
+ if (.not. USE_EXTERNAL_MESH) then
! boundary parameters
open(unit=27,file=prname(1:len_trim(prname))//'ibelm.bin',status='old',action='read',form='unformatted')
read(27) ibelm_xmin
@@ -770,8 +934,308 @@
endif
endif
+
+ endif ! end of (.not. USE_EXTERNAL_MESH)
+! detecting surface points/elements (based on valence check on NGLL points) for external mesh
+ if (USE_EXTERNAL_MESH) then
+ allocate(valence_external_mesh(NGLOB_AB))
+ allocate(ispec_is_surface_external_mesh(NSPEC_AB))
+ allocate(iglob_is_surface_external_mesh(NGLOB_AB))
+ if (.not. RECVS_CAN_BE_BURIED_EXT_MESH) then
+ valence_external_mesh(:) = 0
+ ispec_is_surface_external_mesh(:) = .false.
+ iglob_is_surface_external_mesh(:) = .false.
+ do ispec = 1, NSPEC_AB
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+ valence_external_mesh(iglob) = valence_external_mesh(iglob) + 1
+
+ enddo
+ enddo
+ enddo
+
+ enddo
+
+ allocate(buffer_send_scalar_i_ext_mesh(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh))
+ allocate(buffer_recv_scalar_i_ext_mesh(max_nibool_interfaces_ext_mesh,ninterfaces_ext_mesh))
+
+ call assemble_MPI_scalar_i_ext_mesh(NPROC,NGLOB_AB,valence_external_mesh, &
+ buffer_send_scalar_i_ext_mesh,buffer_recv_scalar_i_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh &
+ )
+
+ do ispec = 1, NSPEC_AB
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ if ( &
+ (k == 1 .or. k == NGLLZ) .and. (j /= 1 .and. j /= NGLLY) .and. (i /= 1 .and. i /= NGLLX) .or. &
+ (j == 1 .or. j == NGLLY) .and. (k /= 1 .and. k /= NGLLZ) .and. (i /= 1 .and. i /= NGLLX) .or. &
+ (i == 1 .or. i == NGLLX) .and. (k /= 1 .and. k /= NGLLZ) .and. (j /= 1 .and. j /= NGLLY) &
+ ) then
+ iglob = ibool(i,j,k,ispec)
+ if (valence_external_mesh(iglob) == 1) then
+ ispec_is_surface_external_mesh(ispec) = .true.
+
+ if (k == 1 .or. k == NGLLZ) then
+ do jj = 1, NGLLY
+ do ii = 1, NGLLX
+ iglob_is_surface_external_mesh(ibool(ii,jj,k,ispec)) = .true.
+ enddo
+ enddo
+ endif
+ if (j == 1 .or. j == NGLLY) then
+ do kk = 1, NGLLZ
+ do ii = 1, NGLLX
+ iglob_is_surface_external_mesh(ibool(ii,j,kk,ispec)) = .true.
+ enddo
+ enddo
+ endif
+ if (i == 1 .or. i == NGLLX) then
+ do kk = 1, NGLLZ
+ do jj = 1, NGLLY
+ iglob_is_surface_external_mesh(ibool(i,jj,kk,ispec)) = .true.
+ enddo
+ enddo
+ endif
+ endif
+
+ endif
+ enddo
+ enddo
+ enddo
+
+ enddo
+
+ if (EXTERNAL_MESH_MOVIE_SURFACE .or. EXTERNAL_MESH_CREATE_SHAKEMAP) then
+ nfaces_surface_external_mesh = 0
+ do ispec = 1, NSPEC_AB
+ iglob = ibool(2,2,1,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ endif
+ iglob = ibool(2,2,NGLLZ,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ endif
+ iglob = ibool(2,1,2,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ endif
+ iglob = ibool(2,NGLLY,2,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ endif
+ iglob = ibool(1,2,2,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ endif
+ iglob = ibool(NGLLX,2,2,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ endif
+ enddo
+
+ allocate(nfaces_perproc_surface_ext_mesh(NPROC))
+ allocate(faces_surface_offset_ext_mesh(NPROC))
+ if (nfaces_surface_external_mesh == 0) then
+ if (USE_HIGHRES_FOR_MOVIES) then
+ allocate(faces_surface_external_mesh(NGLLX*NGLLY,1))
+ allocate(store_val_x_external_mesh(NGLLX*NGLLY*1))
+ allocate(store_val_y_external_mesh(NGLLX*NGLLY*1))
+ allocate(store_val_z_external_mesh(NGLLX*NGLLY*1))
+ allocate(store_val_ux_external_mesh(NGLLX*NGLLY*1))
+ allocate(store_val_uy_external_mesh(NGLLX*NGLLY*1))
+ allocate(store_val_uz_external_mesh(NGLLX*NGLLY*1))
+ else
+ allocate(faces_surface_external_mesh(NGNOD2D,1))
+ allocate(store_val_x_external_mesh(NGNOD2D*1))
+ allocate(store_val_y_external_mesh(NGNOD2D*1))
+ allocate(store_val_z_external_mesh(NGNOD2D*1))
+ allocate(store_val_ux_external_mesh(NGNOD2D*1))
+ allocate(store_val_uy_external_mesh(NGNOD2D*1))
+ allocate(store_val_uz_external_mesh(NGNOD2D*1))
+ endif
+ else
+ if (USE_HIGHRES_FOR_MOVIES) then
+ allocate(faces_surface_external_mesh(NGLLX*NGLLY,nfaces_surface_external_mesh))
+ allocate(store_val_x_external_mesh(NGLLX*NGLLY*nfaces_surface_external_mesh))
+ allocate(store_val_y_external_mesh(NGLLX*NGLLY*nfaces_surface_external_mesh))
+ allocate(store_val_z_external_mesh(NGLLX*NGLLY*nfaces_surface_external_mesh))
+ allocate(store_val_ux_external_mesh(NGLLX*NGLLY*nfaces_surface_external_mesh))
+ allocate(store_val_uy_external_mesh(NGLLX*NGLLY*nfaces_surface_external_mesh))
+ allocate(store_val_uz_external_mesh(NGLLX*NGLLY*nfaces_surface_external_mesh))
+ else
+ allocate(faces_surface_external_mesh(NGNOD2D,nfaces_surface_external_mesh))
+ allocate(store_val_x_external_mesh(NGNOD2D*nfaces_surface_external_mesh))
+ allocate(store_val_y_external_mesh(NGNOD2D*nfaces_surface_external_mesh))
+ allocate(store_val_z_external_mesh(NGNOD2D*nfaces_surface_external_mesh))
+ allocate(store_val_ux_external_mesh(NGNOD2D*nfaces_surface_external_mesh))
+ allocate(store_val_uy_external_mesh(NGNOD2D*nfaces_surface_external_mesh))
+ allocate(store_val_uz_external_mesh(NGNOD2D*nfaces_surface_external_mesh))
+ endif
+ endif
+ call sum_all_i(nfaces_surface_external_mesh,nfaces_surface_glob_ext_mesh)
+ if (USE_HIGHRES_FOR_MOVIES) then
+ allocate(store_val_x_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_y_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_z_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_ux_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_uy_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_uz_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
+ else
+ allocate(store_val_x_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_y_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_z_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_ux_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_uy_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_uz_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
+ endif
+ call gather_all_i(nfaces_surface_external_mesh,1,nfaces_perproc_surface_ext_mesh,1,NPROC)
+
+ faces_surface_offset_ext_mesh(1) = 0
+ do i = 2, NPROC
+ faces_surface_offset_ext_mesh(i) = sum(nfaces_perproc_surface_ext_mesh(1:i-1))
+ enddo
+ if (USE_HIGHRES_FOR_MOVIES) then
+ faces_surface_offset_ext_mesh(:) = faces_surface_offset_ext_mesh(:)*NGLLX*NGLLY
+ else
+ faces_surface_offset_ext_mesh(:) = faces_surface_offset_ext_mesh(:)*NGNOD2D
+ endif
+
+ nfaces_surface_external_mesh = 0
+ do ispec = 1, NSPEC_AB
+ if (ispec_is_surface_external_mesh(ispec)) then
+ iglob = ibool(2,2,1,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ if (USE_HIGHRES_FOR_MOVIES) then
+ ipoin =0
+ do j = NGLLY, 1, -1
+ do i = 1, NGLLX
+ ipoin = ipoin+1
+ faces_surface_external_mesh(ipoin,nfaces_surface_external_mesh) = ibool(i,j,1,ispec)
+ enddo
+ enddo
+ else
+ faces_surface_external_mesh(1,nfaces_surface_external_mesh) = ibool(1,1,1,ispec)
+ faces_surface_external_mesh(2,nfaces_surface_external_mesh) = ibool(1,NGLLY,1,ispec)
+ faces_surface_external_mesh(3,nfaces_surface_external_mesh) = ibool(NGLLX,NGLLY,1,ispec)
+ faces_surface_external_mesh(4,nfaces_surface_external_mesh) = ibool(NGLLX,1,1,ispec)
+ endif
+ endif
+ iglob = ibool(2,2,NGLLZ,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ if (USE_HIGHRES_FOR_MOVIES) then
+ ipoin =0
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ ipoin = ipoin+1
+ faces_surface_external_mesh(ipoin,nfaces_surface_external_mesh) = ibool(i,j,NGLLZ,ispec)
+ enddo
+ enddo
+ else
+ faces_surface_external_mesh(1,nfaces_surface_external_mesh) = ibool(1,1,NGLLZ,ispec)
+ faces_surface_external_mesh(2,nfaces_surface_external_mesh) = ibool(NGLLX,1,NGLLZ,ispec)
+ faces_surface_external_mesh(3,nfaces_surface_external_mesh) = ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ faces_surface_external_mesh(4,nfaces_surface_external_mesh) = ibool(1,NGLLY,NGLLZ,ispec)
+ endif
+ endif
+ iglob = ibool(2,1,2,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ if (USE_HIGHRES_FOR_MOVIES) then
+ ipoin =0
+ do k = 1, NGLLZ
+ do i = 1, NGLLX
+ ipoin = ipoin+1
+ faces_surface_external_mesh(ipoin,nfaces_surface_external_mesh) = ibool(i,1,k,ispec)
+ enddo
+ enddo
+ else
+ faces_surface_external_mesh(1,nfaces_surface_external_mesh) = ibool(1,1,1,ispec)
+ faces_surface_external_mesh(2,nfaces_surface_external_mesh) = ibool(NGLLX,1,1,ispec)
+ faces_surface_external_mesh(3,nfaces_surface_external_mesh) = ibool(NGLLX,1,NGLLZ,ispec)
+ faces_surface_external_mesh(4,nfaces_surface_external_mesh) = ibool(1,1,NGLLZ,ispec)
+ endif
+ endif
+ iglob = ibool(2,NGLLY,2,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ if (USE_HIGHRES_FOR_MOVIES) then
+ ipoin =0
+ do k = 1, NGLLZ
+ do i = NGLLX, 1, -1
+ ipoin = ipoin+1
+ faces_surface_external_mesh(ipoin,nfaces_surface_external_mesh) = ibool(i,NGLLY,k,ispec)
+ enddo
+ enddo
+ else
+ faces_surface_external_mesh(1,nfaces_surface_external_mesh) = ibool(NGLLX,NGLLY,1,ispec)
+ faces_surface_external_mesh(2,nfaces_surface_external_mesh) = ibool(1,NGLLY,1,ispec)
+ faces_surface_external_mesh(3,nfaces_surface_external_mesh) = ibool(1,NGLLY,NGLLZ,ispec)
+ faces_surface_external_mesh(4,nfaces_surface_external_mesh) = ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ endif
+ endif
+ iglob = ibool(1,2,2,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ if (USE_HIGHRES_FOR_MOVIES) then
+ ipoin =0
+ do k = 1, NGLLZ
+ do j = NGLLY, 1, -1
+ ipoin = ipoin+1
+ faces_surface_external_mesh(ipoin,nfaces_surface_external_mesh) = ibool(1,j,k,ispec)
+ enddo
+ enddo
+ else
+ faces_surface_external_mesh(1,nfaces_surface_external_mesh) = ibool(1,NGLLY,1,ispec)
+ faces_surface_external_mesh(2,nfaces_surface_external_mesh) = ibool(1,1,1,ispec)
+ faces_surface_external_mesh(3,nfaces_surface_external_mesh) = ibool(1,1,NGLLZ,ispec)
+ faces_surface_external_mesh(4,nfaces_surface_external_mesh) = ibool(1,NGLLY,NGLLZ,ispec)
+ endif
+ endif
+ iglob = ibool(NGLLX,2,2,ispec)
+ if (iglob_is_surface_external_mesh(iglob)) then
+ nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
+ if (USE_HIGHRES_FOR_MOVIES) then
+ ipoin =0
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ ipoin = ipoin+1
+ faces_surface_external_mesh(ipoin,nfaces_surface_external_mesh) = ibool(NGLLX,j,k,ispec)
+ enddo
+ enddo
+ else
+ faces_surface_external_mesh(1,nfaces_surface_external_mesh) = ibool(NGLLX,1,1,ispec)
+ faces_surface_external_mesh(2,nfaces_surface_external_mesh) = ibool(NGLLX,NGLLY,1,ispec)
+ faces_surface_external_mesh(3,nfaces_surface_external_mesh) = ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ faces_surface_external_mesh(4,nfaces_surface_external_mesh) = ibool(NGLLX,1,NGLLZ,ispec)
+ endif
+ endif
+
+ endif
+ enddo
+
+ if (myrank == 0) then
+ print *, nfaces_perproc_surface_ext_mesh
+ print *, nfaces_surface_glob_ext_mesh
+
+ endif
+
+ endif
+
+ endif !
+
+ endif
+
! $$$$$$$$$$$$$$$$$$$$$$$$ SOURCES $$$$$$$$$$$$$$$$$
! read topography and bathymetry file
@@ -795,6 +1259,11 @@
write(IMAIN,*)
endif
+ else
+ NX_TOPO = 1
+ NY_TOPO = 1
+ allocate(itopo_bathy(NX_TOPO,NY_TOPO))
+
endif
! write source and receiver VTK files for Paraview
@@ -825,6 +1294,7 @@
allocate(hdur_gaussian(NSOURCES))
allocate(utm_x_source(NSOURCES))
allocate(utm_y_source(NSOURCES))
+ allocate(nu_source(3,3,NSOURCES))
! locate sources in the mesh
call locate_source(ibool,NSOURCES,myrank,NSPEC_AB,NGLOB_AB, &
@@ -836,7 +1306,9 @@
LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX,Z_DEPTH_BLOCK, &
TOPOGRAPHY,itopo_bathy,UTM_PROJECTION_ZONE, &
PRINT_SOURCE_TIME_FUNCTION,SUPPRESS_UTM_PROJECTION, &
- NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO)
+ NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO, &
+ nu_source,iglob_is_surface_external_mesh,ispec_is_surface_external_mesh &
+ )
if(minval(t_cmt) /= 0.) call exit_MPI(myrank,'one t_cmt must be zero, others must be positive')
@@ -909,7 +1381,9 @@
xi_receiver,eta_receiver,gamma_receiver,station_name,network_name,nu, &
NPROC,utm_x_source(1),utm_y_source(1), &
TOPOGRAPHY,itopo_bathy,UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION, &
- NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO)
+ NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO, &
+ iglob_is_surface_external_mesh,ispec_is_surface_external_mesh &
+)
!###################### SOURCE ARRAYS ################
@@ -1105,10 +1579,19 @@
call sync_all()
! the mass matrix needs to be assembled with MPI here once and for all
+ if (USE_EXTERNAL_MESH) then
+ call assemble_MPI_scalar_ext_mesh(NPROC,NGLOB_AB,rmass, &
+ buffer_send_scalar_ext_mesh,buffer_recv_scalar_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh &
+ )
+ else
call assemble_MPI_scalar(rmass,iproc_xi,iproc_eta,addressing, &
iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
buffer_send_faces_scalar,buffer_received_faces_scalar,npoin2D_xi,npoin2D_eta, &
NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
+ endif
if(myrank == 0) write(IMAIN,*) 'end assembling MPI mass matrix'
@@ -1427,6 +1910,7 @@
! ************* MAIN LOOP OVER THE TIME STEPS *************
! *********************************************************
+
do it = 1,NSTEP
! compute the maximum of the norm of the displacement
@@ -1507,6 +1991,8 @@
ispec2D_moho_bot = 0
endif
+ if (.not. USE_EXTERNAL_MESH) then
+
do ispec = 1,NSPEC_AB
if (SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3) then
@@ -2404,11 +2890,44 @@
if (SIMULATION_TYPE == 1) then
- do isource = 1,NSOURCES
+!!!!!!!!!! do isource = 1,NSOURCES
+ do isource = 1,1
+ if(FASTER_SOURCES_POINTS_ONLY) then
+
! add the source (only if this proc carries the source)
if(myrank == islice_selected_source(isource)) then
+ iglob = ibool(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec_selected_source(isource))
+ f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
+ t0 = 1.2d0/f0
+if (it == 1 .and. myrank == 0) then
+
+print *,'using a source of dominant frequency ',f0
+print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+endif
+
+
+ ! we use nu_source(:,3) here because we want a source normal to the surface.
+ ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
+ !accel(:,iglob) = accel(:,iglob) + &
+ ! sngl(nu_source(:,3,isource) * 10000000.d0 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
+ ! exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
+accel(:,iglob) = accel(:,iglob) + &
+ sngl(nu_source(:,3,isource) * 1.d10 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
+ exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
+
+ endif
+
+ else
+
+! add the source (only if this proc carries the source)
+ if(myrank == islice_selected_source(isource)) then
+
stf = comp_source_time_function(dble(it-1)*DT-t0-t_cmt(isource),hdur_gaussian(isource))
! distinguish between single and double precision for reals
@@ -2430,6 +2949,8 @@
endif
+ endif ! end of if(FASTER_SOURCES_POINTS_ONLY)
+
enddo
endif
@@ -2490,12 +3011,41 @@
endif
+ endif ! if (.not. USE_EXTERNAL_MESH)
! assemble all the contributions between slices using MPI
+ if (USE_EXTERNAL_MESH) then
+ call compute_forces(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ kappastore,mustore,jacobian,ibool,ispec_is_inner_ext_mesh,.false., &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt &
+ )
+
+ call assemble_MPI_vector_ext_mesh_s(NPROC,NGLOB_AB,accel, &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh &
+ )
+
+ call compute_forces(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ kappastore,mustore,jacobian,ibool,ispec_is_inner_ext_mesh,.true., &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt &
+ )
+
+ call assemble_MPI_vector_ext_mesh_w(NPROC,NGLOB_AB,accel, &
+ buffer_send_vector_ext_mesh,buffer_recv_vector_ext_mesh, &
+ ninterfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
+ request_send_vector_ext_mesh,request_recv_vector_ext_mesh &
+ )
+ else
call assemble_MPI_vector(accel,iproc_xi,iproc_eta,addressing, &
iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
buffer_send_faces_vector,buffer_received_faces_vector,npoin2D_xi,npoin2D_eta, &
NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
+ endif
if (SIMULATION_TYPE == 3) call assemble_MPI_vector(b_accel,iproc_xi,iproc_eta,addressing, &
iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
buffer_send_faces_vector,buffer_received_faces_vector,npoin2D_xi,npoin2D_eta, &
@@ -2591,6 +3141,22 @@
irec = number_receiver_global(irec_local)
! perform the general interpolation using Lagrange polynomials
+ if(FASTER_RECEIVERS_POINTS_ONLY) then
+
+ iglob = ibool(nint(xi_receiver(irec)),nint(eta_receiver(irec)), &
+ nint(gamma_receiver(irec)),ispec_selected_rec(irec))
+ dxd = dble(displ(1,iglob))
+ dyd = dble(displ(2,iglob))
+ dzd = dble(displ(3,iglob))
+ vxd = dble(veloc(1,iglob))
+ vyd = dble(veloc(2,iglob))
+ vzd = dble(veloc(3,iglob))
+ axd = dble(accel(1,iglob))
+ ayd = dble(accel(2,iglob))
+ azd = dble(accel(3,iglob))
+
+ else
+
dxd = ZERO
dyd = ZERO
dzd = ZERO
@@ -2633,6 +3199,8 @@
enddo
enddo
+
+
else if (SIMULATION_TYPE == 2) then
do k = 1,NGLLZ
@@ -2702,8 +3270,8 @@
enddo
endif
+ endif ! end of if(FASTER_RECEIVERS_POINTS_ONLY)
-
! store North, East and Vertical components
! distinguish between single and double precision for reals
@@ -2822,6 +3390,269 @@
endif
endif
+
+ if (EXTERNAL_MESH_CREATE_SHAKEMAP) then
+ if (it == 1) then
+
+ store_val_ux_external_mesh(:) = -HUGEVAL
+ store_val_uy_external_mesh(:) = -HUGEVAL
+ store_val_uz_external_mesh(:) = -HUGEVAL
+ do ispec = 1,nfaces_surface_external_mesh
+ if (USE_HIGHRES_FOR_MOVIES) then
+ do ipoin = 1, NGLLX*NGLLY
+ store_val_x_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = xstore(faces_surface_external_mesh(ipoin,ispec))
+ store_val_y_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = ystore(faces_surface_external_mesh(ipoin,ispec))
+ store_val_z_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = zstore(faces_surface_external_mesh(ipoin,ispec))
+ enddo
+ else
+ store_val_x_external_mesh(NGNOD2D*(ispec-1)+1) = xstore(faces_surface_external_mesh(1,ispec))
+ store_val_x_external_mesh(NGNOD2D*(ispec-1)+2) = xstore(faces_surface_external_mesh(2,ispec))
+ store_val_x_external_mesh(NGNOD2D*(ispec-1)+3) = xstore(faces_surface_external_mesh(3,ispec))
+ store_val_x_external_mesh(NGNOD2D*(ispec-1)+4) = xstore(faces_surface_external_mesh(4,ispec))
+ store_val_y_external_mesh(NGNOD2D*(ispec-1)+1) = ystore(faces_surface_external_mesh(1,ispec))
+ store_val_y_external_mesh(NGNOD2D*(ispec-1)+2) = ystore(faces_surface_external_mesh(2,ispec))
+ store_val_y_external_mesh(NGNOD2D*(ispec-1)+3) = ystore(faces_surface_external_mesh(3,ispec))
+ store_val_y_external_mesh(NGNOD2D*(ispec-1)+4) = ystore(faces_surface_external_mesh(4,ispec))
+ store_val_z_external_mesh(NGNOD2D*(ispec-1)+1) = zstore(faces_surface_external_mesh(1,ispec))
+ store_val_z_external_mesh(NGNOD2D*(ispec-1)+2) = zstore(faces_surface_external_mesh(2,ispec))
+ store_val_z_external_mesh(NGNOD2D*(ispec-1)+3) = zstore(faces_surface_external_mesh(3,ispec))
+ store_val_z_external_mesh(NGNOD2D*(ispec-1)+4) = zstore(faces_surface_external_mesh(4,ispec))
+ endif
+ enddo
+ endif
+
+ do ispec = 1,nfaces_surface_external_mesh
+ if (USE_HIGHRES_FOR_MOVIES) then
+ do ipoin = 1, NGLLX*NGLLY
+ store_val_ux_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = &
+ max(store_val_ux_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin), &
+ sqrt(displ(1,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ displ(2,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ displ(3,faces_surface_external_mesh(ipoin,ispec))**2))
+ store_val_uy_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = &
+ max(store_val_uy_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin), &
+ sqrt(veloc(1,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ veloc(2,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ veloc(3,faces_surface_external_mesh(ipoin,ispec))**2))
+ store_val_uz_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = &
+ max(store_val_uz_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin), &
+ sqrt(accel(1,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ accel(2,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ accel(3,faces_surface_external_mesh(ipoin,ispec))**2))
+
+ enddo
+ else
+ store_val_ux_external_mesh(NGNOD2D*(ispec-1)+1) = &
+ max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+1), &
+ sqrt(displ(1,faces_surface_external_mesh(1,ispec))**2 + &
+ displ(2,faces_surface_external_mesh(1,ispec))**2 + &
+ displ(3,faces_surface_external_mesh(1,ispec))**2))
+ store_val_ux_external_mesh(NGNOD2D*(ispec-1)+1) = &
+ max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+2), &
+ sqrt(displ(1,faces_surface_external_mesh(2,ispec))**2 + &
+ displ(2,faces_surface_external_mesh(2,ispec))**2 + &
+ displ(3,faces_surface_external_mesh(2,ispec))**2))
+ store_val_ux_external_mesh(NGNOD2D*(ispec-1)+3) = &
+ max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+3), &
+ sqrt(displ(1,faces_surface_external_mesh(3,ispec))**2 + &
+ displ(2,faces_surface_external_mesh(3,ispec))**2 + &
+ displ(3,faces_surface_external_mesh(3,ispec))**2))
+ store_val_ux_external_mesh(NGNOD2D*(ispec-1)+4) = &
+ max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+4), &
+ sqrt(displ(1,faces_surface_external_mesh(4,ispec))**2 + &
+ displ(2,faces_surface_external_mesh(4,ispec))**2 + &
+ displ(3,faces_surface_external_mesh(4,ispec))**2))
+ store_val_uy_external_mesh(NGNOD2D*(ispec-1)+1) = &
+ max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+1), &
+ sqrt(veloc(1,faces_surface_external_mesh(1,ispec))**2 + &
+ veloc(2,faces_surface_external_mesh(1,ispec))**2 + &
+ veloc(3,faces_surface_external_mesh(1,ispec))**2))
+ store_val_uy_external_mesh(NGNOD2D*(ispec-1)+1) = &
+ max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+2), &
+ sqrt(veloc(1,faces_surface_external_mesh(2,ispec))**2 + &
+ veloc(2,faces_surface_external_mesh(2,ispec))**2 + &
+ veloc(3,faces_surface_external_mesh(2,ispec))**2))
+ store_val_uy_external_mesh(NGNOD2D*(ispec-1)+3) = &
+ max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+3), &
+ sqrt(veloc(1,faces_surface_external_mesh(3,ispec))**2 + &
+ veloc(2,faces_surface_external_mesh(3,ispec))**2 + &
+ veloc(3,faces_surface_external_mesh(3,ispec))**2))
+ store_val_uy_external_mesh(NGNOD2D*(ispec-1)+4) = &
+ max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+4), &
+ sqrt(veloc(1,faces_surface_external_mesh(4,ispec))**2 + &
+ veloc(2,faces_surface_external_mesh(4,ispec))**2 + &
+ veloc(3,faces_surface_external_mesh(4,ispec))**2))
+ store_val_uz_external_mesh(NGNOD2D*(ispec-1)+1) = &
+ max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+1), &
+ sqrt(accel(1,faces_surface_external_mesh(1,ispec))**2 + &
+ accel(2,faces_surface_external_mesh(1,ispec))**2 + &
+ accel(3,faces_surface_external_mesh(1,ispec))**2))
+ store_val_uz_external_mesh(NGNOD2D*(ispec-1)+1) = &
+ max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+2), &
+ sqrt(accel(1,faces_surface_external_mesh(2,ispec))**2 + &
+ accel(2,faces_surface_external_mesh(2,ispec))**2 + &
+ accel(3,faces_surface_external_mesh(2,ispec))**2))
+ store_val_uz_external_mesh(NGNOD2D*(ispec-1)+3) = &
+ max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+3), &
+ sqrt(accel(1,faces_surface_external_mesh(3,ispec))**2 + &
+ accel(2,faces_surface_external_mesh(3,ispec))**2 + &
+ accel(3,faces_surface_external_mesh(3,ispec))**2))
+ store_val_uz_external_mesh(NGNOD2D*(ispec-1)+4) = &
+ max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+4), &
+ sqrt(accel(1,faces_surface_external_mesh(4,ispec))**2 + &
+ accel(2,faces_surface_external_mesh(4,ispec))**2 + &
+ accel(3,faces_surface_external_mesh(4,ispec))**2))
+ endif
+ enddo
+
+ if (it == NSTEP) then
+ if (USE_HIGHRES_FOR_MOVIES) then
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ else
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ endif
+
+ if(myrank == 0) then
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/shakingdata',status='unknown',form='unformatted')
+ write(IOUT) store_val_x_all_external_mesh
+ write(IOUT) store_val_y_all_external_mesh
+ write(IOUT) store_val_z_all_external_mesh
+ write(IOUT) store_val_ux_all_external_mesh
+ write(IOUT) store_val_uy_all_external_mesh
+ write(IOUT) store_val_uz_all_external_mesh
+ close(IOUT)
+ endif
+ endif
+
+ endif
+
+ if(USE_EXTERNAL_MESH .and. EXTERNAL_MESH_MOVIE_SURFACE .and. mod(it,NTSTEP_BETWEEN_FRAMES) == 0) then
+! get coordinates of surface mesh and surface displacement
+ do ispec = 1,nfaces_surface_external_mesh
+ if (USE_HIGHRES_FOR_MOVIES) then
+ do ipoin = 1, NGLLX*NGLLY
+ store_val_x_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = xstore(faces_surface_external_mesh(ipoin,ispec))
+ store_val_y_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = ystore(faces_surface_external_mesh(ipoin,ispec))
+ store_val_z_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = zstore(faces_surface_external_mesh(ipoin,ispec))
+ store_val_ux_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = veloc(1,faces_surface_external_mesh(ipoin,ispec))
+ store_val_uy_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = veloc(2,faces_surface_external_mesh(ipoin,ispec))
+ store_val_uz_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = veloc(3,faces_surface_external_mesh(ipoin,ispec))
+ enddo
+ else
+ store_val_x_external_mesh(NGNOD2D*(ispec-1)+1) = xstore(faces_surface_external_mesh(1,ispec))
+ store_val_x_external_mesh(NGNOD2D*(ispec-1)+2) = xstore(faces_surface_external_mesh(2,ispec))
+ store_val_x_external_mesh(NGNOD2D*(ispec-1)+3) = xstore(faces_surface_external_mesh(3,ispec))
+ store_val_x_external_mesh(NGNOD2D*(ispec-1)+4) = xstore(faces_surface_external_mesh(4,ispec))
+ store_val_y_external_mesh(NGNOD2D*(ispec-1)+1) = ystore(faces_surface_external_mesh(1,ispec))
+ store_val_y_external_mesh(NGNOD2D*(ispec-1)+2) = ystore(faces_surface_external_mesh(2,ispec))
+ store_val_y_external_mesh(NGNOD2D*(ispec-1)+3) = ystore(faces_surface_external_mesh(3,ispec))
+ store_val_y_external_mesh(NGNOD2D*(ispec-1)+4) = ystore(faces_surface_external_mesh(4,ispec))
+ store_val_z_external_mesh(NGNOD2D*(ispec-1)+1) = zstore(faces_surface_external_mesh(1,ispec))
+ store_val_z_external_mesh(NGNOD2D*(ispec-1)+2) = zstore(faces_surface_external_mesh(2,ispec))
+ store_val_z_external_mesh(NGNOD2D*(ispec-1)+3) = zstore(faces_surface_external_mesh(3,ispec))
+ store_val_z_external_mesh(NGNOD2D*(ispec-1)+4) = zstore(faces_surface_external_mesh(4,ispec))
+ store_val_ux_external_mesh(NGNOD2D*(ispec-1)+1) = veloc(1,faces_surface_external_mesh(1,ispec))
+ store_val_ux_external_mesh(NGNOD2D*(ispec-1)+2) = veloc(1,faces_surface_external_mesh(2,ispec))
+ store_val_ux_external_mesh(NGNOD2D*(ispec-1)+3) = veloc(1,faces_surface_external_mesh(3,ispec))
+ store_val_ux_external_mesh(NGNOD2D*(ispec-1)+4) = veloc(1,faces_surface_external_mesh(4,ispec))
+ store_val_uy_external_mesh(NGNOD2D*(ispec-1)+1) = veloc(2,faces_surface_external_mesh(1,ispec))
+ store_val_uy_external_mesh(NGNOD2D*(ispec-1)+2) = veloc(2,faces_surface_external_mesh(2,ispec))
+ store_val_uy_external_mesh(NGNOD2D*(ispec-1)+3) = veloc(2,faces_surface_external_mesh(3,ispec))
+ store_val_uy_external_mesh(NGNOD2D*(ispec-1)+4) = veloc(2,faces_surface_external_mesh(4,ispec))
+ store_val_uz_external_mesh(NGNOD2D*(ispec-1)+1) = veloc(3,faces_surface_external_mesh(1,ispec))
+ store_val_uz_external_mesh(NGNOD2D*(ispec-1)+2) = veloc(3,faces_surface_external_mesh(2,ispec))
+ store_val_uz_external_mesh(NGNOD2D*(ispec-1)+3) = veloc(3,faces_surface_external_mesh(3,ispec))
+ store_val_uz_external_mesh(NGNOD2D*(ispec-1)+4) = veloc(3,faces_surface_external_mesh(4,ispec))
+ endif
+ enddo
+
+ if (USE_HIGHRES_FOR_MOVIES) then
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
+ store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ else
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_external_mesh*NGNOD2D,&
+ store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ endif
+
+ if(myrank == 0) then
+ write(outputname,"('/moviedata',i6.6)") it
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown',form='unformatted')
+ write(IOUT) store_val_x_all_external_mesh
+ write(IOUT) store_val_y_all_external_mesh
+ write(IOUT) store_val_z_all_external_mesh
+ write(IOUT) store_val_ux_all_external_mesh
+ write(IOUT) store_val_uy_all_external_mesh
+ write(IOUT) store_val_uz_all_external_mesh
+ close(IOUT)
+ endif
+ endif
+
! save MOVIE on the SURFACE
if(MOVIE_SURFACE .and. mod(it,NTSTEP_BETWEEN_FRAMES) == 0) then
@@ -2872,7 +3703,7 @@
enddo
enddo ! ispec_top
endif
-
+
ispec = nmovie_points
call gather_all_cr(store_val_x,ispec,store_val_x_all,ispec,NPROC)
More information about the CIG-COMMITS
mailing list