[cig-commits] r8593 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:59:20 PST 2007
Author: walter
Date: 2007-12-07 15:59:18 -0800 (Fri, 07 Dec 2007)
New Revision: 8593
Modified:
seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90
seismo/2D/SPECFEM2D/trunk/attenuation_model.f90
seismo/2D/SPECFEM2D/trunk/checkgrid.F90
seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90
seismo/2D/SPECFEM2D/trunk/compute_energy.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/compute_gradient_attenuation.f90
seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90
seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90
seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90
seismo/2D/SPECFEM2D/trunk/create_color_image.f90
seismo/2D/SPECFEM2D/trunk/createnum_fast.f90
seismo/2D/SPECFEM2D/trunk/createnum_slow.f90
seismo/2D/SPECFEM2D/trunk/datim.f90
seismo/2D/SPECFEM2D/trunk/define_derivation_matrices.f90
seismo/2D/SPECFEM2D/trunk/define_external_model.f90
seismo/2D/SPECFEM2D/trunk/define_shape_functions.f90
seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90
seismo/2D/SPECFEM2D/trunk/gmat01.f90
seismo/2D/SPECFEM2D/trunk/lagrange_poly.f90
seismo/2D/SPECFEM2D/trunk/locate_receivers.F90
seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90
seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
seismo/2D/SPECFEM2D/trunk/plotgll.f90
seismo/2D/SPECFEM2D/trunk/plotpost.F90
seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90
seismo/2D/SPECFEM2D/trunk/recompute_jacobian.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
Log:
updated the copyright information
Modified: seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/assemble_MPI.F90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -1,5 +1,18 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.2
+! ------------------------------
+!
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
+! University of Pau, France
+!
+! (c) November 2007
+!
+!========================================================================
+
subroutine prepare_assemble_MPI (nspec,ibool, &
- knods, ngnod, &
+ knods, ngnod, &
npoin, elastic, &
ninterface, max_interface_size, &
my_nelmnts_neighbours, my_interfaces, &
@@ -10,10 +23,9 @@
mask_ispec_inner_outer &
)
-
implicit none
+
include 'constants.h'
-
integer, intent(in) :: nspec, npoin, ngnod
logical, dimension(nspec), intent(in) :: elastic
@@ -25,18 +37,17 @@
integer, dimension(ninterface) :: my_nelmnts_neighbours
integer, dimension(4,max_interface_size,ninterface) :: my_interfaces
integer, dimension(NGLLX*max_interface_size,ninterface) :: &
- ibool_interfaces_acoustic,ibool_interfaces_elastic
+ ibool_interfaces_acoustic,ibool_interfaces_elastic
integer, dimension(ninterface) :: &
nibool_interfaces_acoustic,nibool_interfaces_elastic
-
integer, dimension(ninterface), intent(out) :: &
inum_interfaces_acoustic, inum_interfaces_elastic
integer, intent(out) :: ninterface_acoustic, ninterface_elastic
integer :: num_interface
integer :: ispec_interface
-
+
logical, dimension(nspec), intent(inout) :: mask_ispec_inner_outer
logical, dimension(npoin) :: mask_ibool_acoustic
@@ -54,58 +65,57 @@
integer :: npoin_interface_elastic
integer :: ix,iz
-
+
integer :: sens
-
ibool_interfaces_acoustic(:,:) = 0
nibool_interfaces_acoustic(:) = 0
ibool_interfaces_elastic(:,:) = 0
nibool_interfaces_elastic(:) = 0
-
+
do num_interface = 1, ninterface
npoin_interface_acoustic = 0
npoin_interface_elastic = 0
mask_ibool_acoustic(:) = .false.
mask_ibool_elastic(:) = .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, ngnod
- n(k) = knods(k,ispec)
+ n(k) = knods(k,ispec)
end do
e1 = my_interfaces(3,ispec_interface,num_interface)
e2 = my_interfaces(4,ispec_interface,num_interface)
call get_edge(ngnod, n, type, e1, e2, ixmin, ixmax, izmin, izmax, sens)
-
+
do iz = izmin, izmax, sens
do ix = ixmin, ixmax, sens
-
+
if ( elastic(ispec) ) then
-
+
if(.not. mask_ibool_elastic(ibool(ix,iz,ispec))) then
mask_ibool_elastic(ibool(ix,iz,ispec)) = .true.
- npoin_interface_elastic = npoin_interface_elastic + 1
+ npoin_interface_elastic = npoin_interface_elastic + 1
ibool_interfaces_elastic(npoin_interface_elastic,num_interface)=&
ibool(ix,iz,ispec)
end if
else
if(.not. mask_ibool_acoustic(ibool(ix,iz,ispec))) then
mask_ibool_acoustic(ibool(ix,iz,ispec)) = .true.
- npoin_interface_acoustic = npoin_interface_acoustic + 1
+ npoin_interface_acoustic = npoin_interface_acoustic + 1
ibool_interfaces_acoustic(npoin_interface_acoustic,num_interface)=&
ibool(ix,iz,ispec)
end if
end if
end do
end do
-
+
end do
nibool_interfaces_acoustic(num_interface) = npoin_interface_acoustic
nibool_interfaces_elastic(num_interface) = npoin_interface_elastic
-
+
do ispec = 1, nspec
do iz = 1, NGLLZ
do ix = 1, NGLLX
@@ -118,9 +128,7 @@
enddo
enddo
-
end do
-
ninterface_acoustic = 0
ninterface_elastic = 0
@@ -138,10 +146,10 @@
end subroutine prepare_assemble_MPI
-
subroutine get_edge ( ngnod, n, type, e1, e2, ixmin, ixmax, izmin, izmax, sens )
implicit none
+
include "constants.h"
integer, intent(in) :: ngnod
@@ -180,7 +188,7 @@
if ( e1 == n(1) ) then
ixmin = 1
izmin = 1
- if ( e2 == n(2) ) then
+ if ( e2 == n(2) ) then
ixmax = NGLLX
izmax = 1
sens = 1
@@ -234,12 +242,10 @@
end if
end if
end if
-
end subroutine get_edge
-
#ifdef USE_MPI
subroutine create_MPI_req_SEND_RECV_ac( &
@@ -258,7 +264,6 @@
include 'constants.h'
include 'mpif.h'
include 'precision_mpi.h'
-
integer, intent(in) :: ninterface, ninterface_acoustic
integer, dimension(ninterface), intent(in) :: inum_interfaces_acoustic
@@ -271,15 +276,13 @@
integer, dimension(ninterface), intent(in) :: nibool_interfaces_acoustic
integer, dimension(ninterface), intent(in) :: my_neighbours
-
integer :: inum_interface,num_interface
integer :: ier
-
do inum_interface = 1, ninterface_acoustic
-
+
num_interface = inum_interfaces_acoustic(inum_interface)
-
+
call MPI_Send_init ( buffer_send_faces_vector_ac(1,inum_interface), &
nibool_interfaces_acoustic(num_interface), CUSTOM_MPI_TYPE, &
my_neighbours(num_interface), 12, MPI_COMM_WORLD, &
@@ -290,11 +293,9 @@
tab_requests_send_recv_acoustic(ninterface_acoustic+inum_interface), ier)
end do
-
end subroutine create_MPI_req_SEND_RECV_ac
-
subroutine create_MPI_req_SEND_RECV_el( &
ninterface, ninterface_elastic, &
nibool_interfaces_elastic, &
@@ -310,7 +311,7 @@
include 'constants.h'
include 'mpif.h'
- include 'precision_mpi.h'
+ include 'precision_mpi.h'
integer, intent(in) :: ninterface, ninterface_elastic
@@ -326,13 +327,11 @@
integer :: inum_interface,num_interface
integer :: ier
-
-
do inum_interface = 1, ninterface_elastic
-
+
num_interface = inum_interfaces_elastic(inum_interface)
-
+
call MPI_Send_init ( buffer_send_faces_vector_el(1,inum_interface), &
NDIM*nibool_interfaces_elastic(num_interface), CUSTOM_MPI_TYPE, &
my_neighbours(num_interface), 13, MPI_COMM_WORLD, &
@@ -343,11 +342,9 @@
tab_requests_send_recv_elastic(ninterface_elastic+inum_interface), ier)
end do
-
end subroutine create_MPI_req_SEND_RECV_el
-
subroutine assemble_MPI_scalar(myrank,array_val1, array_val2,npoin, &
ninterface, max_interface_size, max_ibool_interfaces_size_ac, max_ibool_interfaces_size_el, &
ibool_interfaces_acoustic,ibool_interfaces_elastic, nibool_interfaces_acoustic,nibool_interfaces_elastic, my_neighbours)
@@ -357,7 +354,6 @@
include 'constants.h'
include 'mpif.h'
-
! array to assemble
real(kind=CUSTOM_REAL), dimension(npoin), intent(inout) :: array_val1, array_val2
@@ -371,7 +367,6 @@
integer, dimension(ninterface), intent(in) :: nibool_interfaces_acoustic,nibool_interfaces_elastic
integer, dimension(ninterface), intent(in) :: my_neighbours
-
double precision, dimension(max_ibool_interfaces_size_ac+max_ibool_interfaces_size_el, ninterface) :: &
buffer_send_faces_scalar, &
buffer_recv_faces_scalar
@@ -379,32 +374,30 @@
integer, dimension(ninterface) :: msg_requests
integer :: ipoin, num_interface
integer :: ier
-
+
integer :: i
-
do num_interface = 1, ninterface
-
+
ipoin = 0
do i = 1, nibool_interfaces_acoustic(num_interface)
ipoin = ipoin + 1
buffer_send_faces_scalar(ipoin,num_interface) = &
array_val1(ibool_interfaces_acoustic(i,num_interface))
end do
-
- do i = 1, nibool_interfaces_elastic(num_interface)
+
+ do i = 1, nibool_interfaces_elastic(num_interface)
ipoin = ipoin + 1
buffer_send_faces_scalar(ipoin,num_interface) = &
array_val2(ibool_interfaces_elastic(i,num_interface))
end do
-
+
call MPI_isend ( buffer_send_faces_scalar(1,num_interface), &
nibool_interfaces_acoustic(num_interface)+nibool_interfaces_elastic(num_interface), MPI_DOUBLE_PRECISION, &
my_neighbours(num_interface), 11, &
MPI_COMM_WORLD, msg_requests(num_interface), ier)
-
+
end do
-
do num_interface = 1, ninterface
call MPI_recv ( buffer_recv_faces_scalar(1,num_interface), &
@@ -418,7 +411,7 @@
array_val1(ibool_interfaces_acoustic(i,num_interface)) = array_val1(ibool_interfaces_acoustic(i,num_interface)) + &
buffer_recv_faces_scalar(ipoin,num_interface)
end do
-
+
do i = 1, nibool_interfaces_elastic(num_interface)
ipoin = ipoin + 1
array_val2(ibool_interfaces_elastic(i,num_interface)) = array_val2(ibool_interfaces_elastic(i,num_interface)) + &
@@ -429,11 +422,9 @@
call MPI_BARRIER(mpi_comm_world,ier)
-
end subroutine assemble_MPI_scalar
-
subroutine assemble_MPI_vector_ac_start(array_val1,npoin, &
ninterface, ninterface_acoustic, &
inum_interfaces_acoustic, &
@@ -448,11 +439,9 @@
include 'constants.h'
include 'mpif.h'
-
! array to assemble
real(kind=CUSTOM_REAL), dimension(npoin), intent(in) :: array_val1
-
integer, intent(in) :: npoin
integer, intent(in) :: ninterface, ninterface_acoustic
integer, dimension(ninterface), intent(in) :: inum_interfaces_acoustic
@@ -466,37 +455,34 @@
integer :: ipoin, num_interface, inum_interface
integer :: ier
-
+
integer :: i
-
do inum_interface = 1, ninterface_acoustic
-
+
num_interface = inum_interfaces_acoustic(inum_interface)
ipoin = 0
do i = 1, nibool_interfaces_acoustic(num_interface)
- ipoin = ipoin + 1
+ ipoin = ipoin + 1
buffer_send_faces_vector_ac(ipoin,inum_interface) = &
array_val1(ibool_interfaces_acoustic(i,num_interface))
end do
-
+
end do
-
+
do inum_interface = 1, ninterface_acoustic*2
call MPI_START(tab_requests_send_recv_acoustic(inum_interface), ier)
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_start unsuccessful in assemble_MPI_vector_start')
end if
end do
-
+
!call MPI_Startall ( ninterface*2, tab_requests_send_recv(1), ier )
-
end subroutine assemble_MPI_vector_ac_start
-
subroutine assemble_MPI_vector_el_start(array_val2,npoin, &
ninterface, ninterface_elastic, &
inum_interfaces_elastic, &
@@ -511,11 +497,9 @@
include 'constants.h'
include 'mpif.h'
-
! array to assemble
real(kind=CUSTOM_REAL), dimension(NDIM,npoin), intent(in) :: array_val2
-
integer, intent(in) :: npoin
integer, intent(in) :: ninterface, ninterface_elastic
integer, dimension(ninterface), intent(in) :: inum_interfaces_elastic
@@ -526,7 +510,6 @@
integer, dimension(ninterface_elastic*2), intent(inout) :: tab_requests_send_recv_elastic
real(CUSTOM_REAL), dimension(max_ibool_interfaces_size_el,ninterface_elastic), intent(inout) :: &
buffer_send_faces_vector_el
-
integer :: ipoin, num_interface, inum_interface
@@ -536,7 +519,7 @@
do inum_interface = 1, ninterface_elastic
-
+
num_interface = inum_interfaces_elastic(inum_interface)
ipoin = 0
@@ -546,23 +529,20 @@
ipoin = ipoin + 2
end do
-
end do
-
+
do inum_interface = 1, ninterface_elastic*2
call MPI_START(tab_requests_send_recv_elastic(inum_interface), ier)
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_start unsuccessful in assemble_MPI_vector_start')
end if
end do
-
+
!call MPI_Startall ( ninterface*2, tab_requests_send_recv(1), ier )
-
end subroutine assemble_MPI_vector_el_start
-
subroutine assemble_MPI_vector_ac_wait(array_val1,npoin, &
ninterface, ninterface_acoustic, &
inum_interfaces_acoustic, &
@@ -577,11 +557,9 @@
include 'constants.h'
include 'mpif.h'
-
! array to assemble
real(kind=CUSTOM_REAL), dimension(npoin), intent(inout) :: array_val1
-
integer, intent(in) :: npoin
integer, intent(in) :: ninterface, ninterface_acoustic
integer, dimension(ninterface), intent(in) :: inum_interfaces_acoustic
@@ -599,14 +577,13 @@
integer :: i
-
call MPI_Waitall ( ninterface_acoustic*2, tab_requests_send_recv_acoustic(1), tab_statuses_acoustic(1,1), ier )
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_WAITALL unsuccessful in assemble_MPI_vector_wait')
end if
-
+
do inum_interface = 1, ninterface_acoustic
-
+
num_interface = inum_interfaces_acoustic(inum_interface)
ipoin = 0
@@ -615,14 +592,12 @@
array_val1(ibool_interfaces_acoustic(i,num_interface)) = array_val1(ibool_interfaces_acoustic(i,num_interface)) + &
buffer_recv_faces_vector_ac(ipoin,inum_interface)
end do
-
+
end do
-
end subroutine assemble_MPI_vector_ac_wait
-
subroutine assemble_MPI_vector_el_wait(array_val2,npoin, &
ninterface, ninterface_elastic, &
inum_interfaces_elastic, &
@@ -637,11 +612,9 @@
include 'constants.h'
include 'mpif.h'
-
! array to assemble
real(kind=CUSTOM_REAL), dimension(NDIM,npoin), intent(inout) :: array_val2
-
integer, intent(in) :: npoin
integer, intent(in) :: ninterface, ninterface_elastic
integer, dimension(ninterface), intent(in) :: inum_interfaces_elastic
@@ -659,14 +632,13 @@
integer :: i
-
call MPI_Waitall ( ninterface_elastic*2, tab_requests_send_recv_elastic(1), tab_statuses_elastic(1,1), ier )
if ( ier /= MPI_SUCCESS ) then
call exit_mpi('MPI_WAITALL unsuccessful in assemble_MPI_vector_wait')
end if
-
+
do inum_interface = 1, ninterface_elastic
-
+
num_interface = inum_interfaces_elastic(inum_interface)
ipoin = 0
@@ -675,16 +647,14 @@
buffer_recv_faces_vector_el(ipoin+1:ipoin+2,inum_interface)
ipoin = ipoin + 2
end do
-
+
end do
-
end subroutine assemble_MPI_vector_el_wait
#endif
-
subroutine exit_MPI(error_msg)
implicit none
@@ -712,15 +682,9 @@
#ifdef USE_MPI
call MPI_ABORT(MPI_COMM_WORLD,30,ier)
call MPI_FINALIZE(ier)
-
+
#endif
stop 'error, program ended in exit_MPI'
-
end subroutine exit_MPI
-
-
-
-
-
Modified: seismo/2D/SPECFEM2D/trunk/attenuation_model.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/attenuation_model.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/attenuation_model.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -30,15 +30,15 @@
double precision, dimension(N_SLS) :: tau_epsilon_nu1,tau_sigma_nu1,tau_epsilon_nu2,tau_sigma_nu2
double precision :: f1_attenuation, f2_attenuation
-
-
+
+
! f1 and f2 are computed as : f2/f1=12 and (log(f1)+log(f2))/2 = log(f0)
f1_attenuation = exp(log(f0_attenuation)-log(12.d0)/2.d0)
f2_attenuation = 12.d0 * f1_attenuation
-! Call of C function that computes attenuation parameters (function in file "attenuation_compute_param.c";
+! Call of C function that computes attenuation parameters (function in file "attenuation_compute_param.c";
! a main can be found in UTILS/attenuation directory).
-! Beware of underscores in this function name; depending on your compiler and compilation options, you will have to add or
+! Beware of underscores in this function name; depending on your compiler and compilation options, you will have to add or
! delete underscores. Also look in file "attenuation_compute_param.c" for this issue.
call attenuation_compute_param(N_SLS, Qp_attenuation, Qs_attenuation, &
f1_attenuation,f2_attenuation, &
Modified: seismo/2D/SPECFEM2D/trunk/checkgrid.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -79,7 +79,7 @@
integer :: num_ispec
integer :: myrank, iproc, nproc
integer :: ier
-
+
#ifdef USE_MPI
integer, dimension(MPI_STATUS_SIZE) :: request_mpi_status
#endif
@@ -1426,7 +1426,7 @@
enddo
any_elastic_glob = any_elastic
-#ifdef USE_MPI
+#ifdef USE_MPI
call MPI_ALLREDUCE (vpmin, vpmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
call MPI_ALLREDUCE (vpmax, vpmax_glob, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ier)
call MPI_ALLREDUCE (vsmin, vsmin_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ier)
@@ -1455,7 +1455,7 @@
lambdaPmax = lambdaPmax_glob
lambdaSmin = lambdaSmin_glob
lambdaSmax = lambdaSmax_glob
-
+
#endif
if ( myrank == 0 ) then
@@ -1478,8 +1478,8 @@
write(IOUT,*)
write(IOUT,*) '*** Max stability for P wave velocity = ',courant_stability_number_max
write(IOUT,*)
-
+
! only if time source is not a Dirac or Heaviside (otherwise maximum frequency of spectrum undefined)
! and if source is not an initial field, for the same reason
if(.not. initialfield .and. time_function_type /= 4 .and. time_function_type /= 5) then
@@ -1535,7 +1535,7 @@
xmax = xmax_glob
zmin = zmin_glob
zmax = zmax_glob
-
+
#endif
! ratio of physical page size/size of the domain meshed
@@ -1653,7 +1653,7 @@
end if
do ispec = 1, nspec
-
+
if ( myrank == 0 ) then
num_ispec = num_ispec + 1
write(24,*) '% elem ',ispec
@@ -1813,14 +1813,14 @@
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
allocate(coorg_recv(2,nspec_recv*5))
allocate(RGB_recv(nspec_recv))
call MPI_RECV (coorg_recv(1,1), nspec_recv*5*2, MPI_DOUBLE_PRECISION, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
call MPI_RECV (RGB_recv(1), nspec_recv, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
do ispec = 1, nspec_recv
num_ispec = num_ispec + 1
write(24,*) '% elem ',num_ispec
@@ -1839,9 +1839,9 @@
end do
deallocate(coorg_recv)
deallocate(RGB_recv)
-
+
end do
-
+
else
call MPI_SEND (nspec, 1, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
call MPI_SEND (coorg_send, nspec*5*2, MPI_DOUBLE_PRECISION, 0, 42, MPI_COMM_WORLD, ier)
@@ -1864,7 +1864,7 @@
!
!--------------------------------------------------------------------------------
!
-
+
if ( myrank == 0 ) then
print *
print *,'Creating PostScript file with mesh dispersion'
@@ -1979,7 +1979,7 @@
write(24,*) '% spectral element mesh'
write(24,*) '%'
write(24,*) '0 setgray'
-
+
num_ispec = 0
end if
@@ -2069,8 +2069,8 @@
coorg_send(2,(ispec-1)*5+5) = z2
end if
-
+
material = kmato(ispec)
mu = elastcoef(2,material)
@@ -2149,9 +2149,9 @@
else
! do not color the elements if not close to the threshold
- if ( myrank == 0 ) then
+ if ( myrank == 0 ) then
write(24,*) 'ST'
- else
+ else
RGB_send(ispec) = 0
end if
endif
@@ -2182,7 +2182,7 @@
else if(lambdaP_local <= 1.20 * lambdaPmin) then
if ( myrank == 0 ) then
write(24,*) '0 0 1 RG GF 0 setgray ST'
- else
+ else
RGB_send(ispec) = 3
end if
@@ -2201,14 +2201,14 @@
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
allocate(coorg_recv(2,nspec_recv*5))
allocate(RGB_recv(nspec_recv))
call MPI_RECV (coorg_recv(1,1), nspec_recv*5*2, MPI_DOUBLE_PRECISION, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
call MPI_RECV (RGB_recv(1), nspec_recv, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
do ispec = 1, nspec_recv
num_ispec = num_ispec + 1
write(24,*) '% elem ',num_ispec
@@ -2228,13 +2228,13 @@
if ( RGB_recv(ispec) == 0) then
write(24,*) 'ST'
end if
-
+
end do
deallocate(coorg_recv)
deallocate(RGB_recv)
end do
-
+
else
call MPI_SEND (nspec, 1, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
call MPI_SEND (coorg_send, nspec*5*2, MPI_DOUBLE_PRECISION, 0, 42, MPI_COMM_WORLD, ier)
@@ -2248,9 +2248,9 @@
write(24,*) '%'
write(24,*) 'grestore'
write(24,*) 'showpage'
-
+
close(24)
-
+
print *,'End of creation of PostScript file with mesh dispersion'
end if
@@ -2367,7 +2367,7 @@
num_ispec = 0
end if
-
+
do ispec = 1, nspec
if ( myrank == 0 ) then
num_ispec = num_ispec + 1
@@ -2412,7 +2412,7 @@
coorg_send(1,(ispec-1)*5+2) = x2
coorg_send(2,(ispec-1)*5+2) = z2
end if
-
+
ir=pointsdisp
is=pointsdisp
x2 = (xinterp(ir,is)-xmin)*ratio_page + orig_x
@@ -2452,7 +2452,7 @@
coorg_send(1,(ispec-1)*5+5) = x2
coorg_send(2,(ispec-1)*5+5) = z2
end if
-
+
if((vpmax-vpmin)/vpmin > 0.02d0) then
if(assign_external_model) then
! use lower-left corner
@@ -2479,21 +2479,21 @@
! display P-velocity model using gray levels
if ( myrank == 0 ) then
write(24,*) sngl(x1),' setgray GF 0 setgray ST'
- else
+ else
greyscale_send(ispec) = sngl(x1)
end if
enddo ! end of loop on all the spectral elements
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
allocate(coorg_recv(2,nspec_recv*5))
allocate(greyscale_recv(nspec_recv))
call MPI_RECV (coorg_recv(1,1), nspec_recv*5*2, MPI_DOUBLE_PRECISION, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
call MPI_RECV (greyscale_recv(1), nspec_recv, MPI_REAL, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
do ispec = 1, nspec_recv
num_ispec = num_ispec + 1
write(24,*) '% elem ',num_ispec
@@ -2505,13 +2505,13 @@
write(24,681) coorg_recv(1,(ispec-1)*5+5), coorg_recv(2,(ispec-1)*5+5)
write(24,*) 'CO'
write(24,*) greyscale_recv(ispec), ' setgray GF 0 setgray ST'
-
+
end do
deallocate(coorg_recv)
deallocate(greyscale_recv)
end do
-
+
else
call MPI_SEND (nspec, 1, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
call MPI_SEND (coorg_send, nspec*5*2, MPI_DOUBLE_PRECISION, 0, 42, MPI_COMM_WORLD, ier)
@@ -2524,14 +2524,14 @@
write(24,*) '%'
write(24,*) 'grestore'
write(24,*) 'showpage'
-
+
close(24)
-
+
print *,'End of creation of PostScript file with velocity model'
end if
-
+
if ( myrank == 0 ) then
print *
print *,'Creating PostScript file with partitioning'
@@ -2642,7 +2642,7 @@
end if
do ispec = 1, nspec
-
+
if ( myrank == 0 ) then
num_ispec = num_ispec + 1
write(24,*) '% elem ',ispec
@@ -2728,25 +2728,25 @@
coorg_send(2,(ispec-1)*5+5) = z2
end if
-
+
if ( myrank == 0 ) then
write(24,*) red(1), green(1), blue(1), 'RG GF 0 setgray ST'
end if
-
+
enddo ! end of loop on all the spectral elements
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
! use a different color for each material set
icol = mod(iproc, NUM_COLORS) + 1
-
+
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
allocate(coorg_recv(2,nspec_recv*5))
call MPI_RECV (coorg_recv(1,1), nspec_recv*5*2, MPI_DOUBLE_PRECISION, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
do ispec = 1, nspec_recv
num_ispec = num_ispec + 1
write(24,*) '% elem ',num_ispec
@@ -2757,18 +2757,18 @@
write(24,681) coorg_recv(1,(ispec-1)*5+4), coorg_recv(2,(ispec-1)*5+4)
write(24,681) coorg_recv(1,(ispec-1)*5+5), coorg_recv(2,(ispec-1)*5+5)
write(24,*) 'CO'
-
+
write(24,*) red(icol), green(icol), blue(icol), ' RG GF 0 setgray ST'
-
+
end do
deallocate(coorg_recv)
-
+
end do
-
+
else
call MPI_SEND (nspec, 1, MPI_INTEGER, 0, 42, MPI_COMM_WORLD, ier)
call MPI_SEND (coorg_send, nspec*5*2, MPI_DOUBLE_PRECISION, 0, 42, MPI_COMM_WORLD, ier)
-
+
end if
#endif
Modified: seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/compute_arrays_source.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -28,12 +28,12 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLX) :: vector_field_element
! pressure in an element
- integer :: N_SLS
+ integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: pressure_element
-
+
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
double precision :: Mu_nu1,Mu_nu2
-
+
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_acoustic,potential_dot_dot_acoustic
logical :: TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -145,7 +145,7 @@
enddo ! end of loop over all spectral elements
-! only for the first call to compute_forces_acoustic (during computation on outer elements)
+! only for the first call to compute_forces_acoustic (during computation on outer elements)
if ( num_phase_inner_outer ) then
!
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -106,7 +106,7 @@
! loop over spectral elements
do ispec_inner_outer = 1,nspec_inner_outer
-
+
! get global numbering for inner or outer elements
ispec = ispec_inner_outer_to_glob(ispec_inner_outer)
Modified: seismo/2D/SPECFEM2D/trunk/compute_gradient_attenuation.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_gradient_attenuation.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/compute_gradient_attenuation.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -1,3 +1,16 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.2
+! ------------------------------
+!
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
+! University of Pau, France
+!
+! (c) November 2007
+!
+!========================================================================
+
subroutine construct_acoustic_surface ( nspec, ngnod, knods, nsurface, surface, tab_surface )
implicit none
@@ -33,15 +46,11 @@
tab_surface(4,i) = izmin
tab_surface(5,i) = izmax
-
enddo
-
end subroutine construct_acoustic_surface
-
-
subroutine get_acoustic_edge ( ngnod, n, type, e1, e2, ixmin, ixmax, izmin, izmax )
implicit none
@@ -142,7 +151,5 @@
endif
endif
-
end subroutine get_acoustic_edge
-
Modified: seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/create_color_image.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/create_color_image.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/create_color_image.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/createnum_fast.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/createnum_fast.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/createnum_fast.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/createnum_slow.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/createnum_slow.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/createnum_slow.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/datim.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/datim.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/datim.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/define_derivation_matrices.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/define_derivation_matrices.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/define_derivation_matrices.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/define_external_model.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/define_external_model.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/define_external_model.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/define_shape_functions.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/define_shape_functions.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/define_shape_functions.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/gmat01.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/gmat01.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/lagrange_poly.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/lagrange_poly.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/lagrange_poly.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/locate_receivers.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_receivers.F90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/locate_receivers.F90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -27,7 +27,7 @@
#endif
integer nrec,nspec,npoin,ngnod,npgeo
- integer, intent(in) :: nproc, myrank
+ integer, intent(in) :: nproc, myrank
integer knods(ngnod,nspec)
double precision coorg(NDIM,npgeo)
@@ -64,14 +64,14 @@
double precision, dimension(nrec) :: st_xval,st_zval
-
+
double precision, dimension(nrec,nproc) :: gather_final_distance
double precision, dimension(nrec,nproc) :: gather_xi_receiver, gather_gamma_receiver
integer, dimension(nrec,nproc) :: gather_ispec_selected_rec
integer, dimension(nrec), intent(inout) :: which_proc_receiver
integer :: ierror
-
+
ierror = 0
! **************
@@ -202,7 +202,7 @@
if ( myrank == 0 ) then
do irec = 1, nrec
which_proc_receiver(irec:irec) = minloc(gather_final_distance(irec,:)) - 1
-
+
end do
end if
@@ -259,7 +259,7 @@
write(IOUT,*)
write(IOUT,*) 'end of receiver detection'
write(IOUT,*)
-
+
end if
@@ -268,7 +268,7 @@
#ifdef USE_MPI
call MPI_BARRIER(MPI_COMM_WORLD,ierror)
-#endif
+#endif
end subroutine locate_receivers
Modified: seismo/2D/SPECFEM2D/trunk/locate_source_force.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_force.F90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -98,40 +98,40 @@
#endif
-
+
! check if this process contains the source
- if ( dist_glob == distminmax ) then
+ if ( dist_glob == distminmax ) then
is_proc_source = 1
end if
-
+
#ifdef USE_MPI
! determining the number of processes that contain the source (useful when the source is located on an interface)
call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierror)
-
+
#else
nb_proc_source = is_proc_source
-
-#endif
-
+
+#endif
+
if ( nb_proc_source < 1 ) then
call exit_MPI('error locating force source')
end if
-
- if ( is_proc_source == 1 ) then
+
+ if ( is_proc_source == 1 ) then
write(iout,200)
-
+
write(iout,"(1x,f12.3,1x,f12.3,1x,f12.3,1x,f12.3,f12.3,1x,i5.5)") x_source,z_source, &
coord(1,iglob_source),coord(2,iglob_source),distmin,nb_proc_source
write(iout,*)
write(iout,*)
write(iout,"('Maximum distance between asked and real =',f12.3)") distminmax
-
+
end if
#ifdef USE_MPI
call MPI_BARRIER(MPI_COMM_WORLD,ierror)
-#endif
+#endif
200 format(//1x,48('=')/,' = S o u r c e s ', &
'r e a l p o s i t i o n s ='/1x,48('=')// &
Modified: seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/locate_source_moment_tensor.F90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -67,7 +67,7 @@
write(IOUT,*) ' locating moment-tensor source'
write(IOUT,*) '*******************************'
write(IOUT,*)
- end if
+ end if
! set distance to huge initial value
distmin=HUGEVAL
@@ -75,7 +75,7 @@
is_proc_source = 0
do ispec=1,nspec
-
+
! loop only on points inside the element
! exclude edges to ensure this point is not shared with other elements
do j=2,NGLLZ-1
@@ -83,7 +83,7 @@
iglob = ibool(i,j,ispec)
dist = sqrt((x_source-dble(coord(1,iglob)))**2 + (z_source-dble(coord(2,iglob)))**2)
-
+
! keep this point if it is closer to the source
if(dist < distmin) then
distmin = dist
@@ -91,56 +91,56 @@
ix_initial_guess = i
iz_initial_guess = j
endif
-
+
enddo
enddo
! end of loop on all the spectral elements
enddo
-
+
#ifdef USE_MPI
! global minimum distance computed over all processes
call MPI_ALLREDUCE (distmin, dist_glob, 1, MPI_DOUBLE_PRECISION, MPI_MIN, MPI_COMM_WORLD, ierror)
-
+
#else
dist_glob = distmin
-#endif
+#endif
! check if this process contains the source
- if ( dist_glob == distmin ) then
+ if ( dist_glob == distmin ) then
is_proc_source = 1
end if
-
-
+
+
#ifdef USE_MPI
! determining the number of processes that contain the source (useful when the source is located on an interface)
call MPI_ALLREDUCE (is_proc_source, nb_proc_source, 1, MPI_INTEGER, MPI_SUM, MPI_COMM_WORLD, ierror)
-
+
#else
nb_proc_source = is_proc_source
-
-#endif
-
+#endif
+
+
#ifdef USE_MPI
! when several processes contain the source, we elect one of them (minimum rank).
if ( nb_proc_source > 1 ) then
-
+
call MPI_ALLGATHER(is_proc_source, 1, MPI_INTEGER, allgather_is_proc_source(1), 1, MPI_INTEGER, MPI_COMM_WORLD, ierror)
locate_is_proc_source = maxloc(allgather_is_proc_source) - 1
-
+
if ( myrank /= locate_is_proc_source(1) ) then
is_proc_source = 0
end if
nb_proc_source = 1
-
+
end if
-
-#endif
+#endif
+
! ****************************************
! find the best (xi,gamma) for each source
! ****************************************
@@ -193,9 +193,9 @@
if ( is_proc_source == 1 ) then
write(IOUT,*)
write(IOUT,*) 'Moment-tensor source:'
-
+
if(final_distance == HUGEVAL) call exit_MPI('error locating moment-tensor source')
-
+
write(IOUT,*) ' original x: ',sngl(x_source)
write(IOUT,*) ' original z: ',sngl(z_source)
write(IOUT,*) 'closest estimate found: ',sngl(final_distance),' m away'
@@ -210,7 +210,7 @@
#ifdef USE_MPI
call MPI_BARRIER(MPI_COMM_WORLD,ierror)
-#endif
+#endif
end subroutine locate_source_moment_tensor
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -79,7 +79,7 @@
integer ixdebregion,ixfinregion,izdebregion,izfinregion
integer iregion,imaterial,nbregion,nb_materials
integer NTSTEP_BETWEEN_OUTPUT_INFO,pointsdisp,subsamp,seismotype,imagetype
- logical generate_STATIONS
+ logical generate_STATIONS
integer ngnod,nt,nx,nz,nxread,nzread,icodematread,ireceiverlines,nreceiverlines
integer, dimension(:), allocatable :: nrec
@@ -114,19 +114,19 @@
! parameters for external mesh
logical :: read_external_mesh
character(len=256) :: mesh_file, nodes_coords_file, materials_file, free_surface_file, absorbing_surface_file
-
-! variables used for storing info about the mesh and partitions
+
+! variables used for storing info about the mesh and partitions
integer, dimension(:), pointer :: elmnts
integer, dimension(:), pointer :: elmnts_bis
double precision, dimension(:,:), pointer :: nodes_coords
integer, dimension(:,:), pointer :: acoustic_surface
integer, dimension(:,:), pointer :: abs_surface
-
+
integer, dimension(:), pointer :: xadj
integer, dimension(:), pointer :: adjncy
integer, dimension(:), pointer :: nnodes_elmnts
integer, dimension(:), pointer :: nodes_elmnts
-
+
integer, dimension(:), pointer :: vwgt
integer, dimension(:), pointer :: adjwgt
integer, dimension(:), pointer :: part
@@ -138,7 +138,7 @@
integer, dimension(:), pointer :: tab_size_interfaces, tab_interfaces
integer, dimension(:), allocatable :: my_interfaces
integer, dimension(:), allocatable :: my_nb_interfaces
-
+
integer :: nedges_coupled, nedges_coupled_loc
integer, dimension(:,:), pointer :: edges_coupled
@@ -166,7 +166,7 @@
character(len=256) :: prname
! variables used for attenuation
- integer :: N_SLS
+ integer :: N_SLS
double precision :: Qp_attenuation
double precision :: Qs_attenuation
double precision :: f0_attenuation
@@ -207,22 +207,22 @@
! read info about partitionning
call read_value_integer(IIN,IGNORE_JUNK,nproc)
if ( nproc <= 0 ) then
- print *, 'Number of processes (nproc) must be greater than or equal to one.'
- stop
+ print *, 'Number of processes (nproc) must be greater than or equal to one.'
+ stop
endif
-
+
#ifndef USE_MPI
if ( nproc > 1 ) then
- print *, 'Number of processes (nproc) must be equal to one when not using MPI.'
+ print *, 'Number of processes (nproc) must be equal to one when not using MPI.'
print *, 'Please recompile with -DUSE_MPI in order to enable use of MPI.'
- stop
+ stop
endif
-
+
#endif
-
+
call read_value_integer(IIN,IGNORE_JUNK,partitionning_method)
call read_value_string(IIN,IGNORE_JUNK,partitionning_strategy)
- select case(partitionning_method)
+ select case(partitionning_method)
case(1)
case(2)
partitionning_strategy = trim(partitionning_strategy)
@@ -233,17 +233,17 @@
metis_options = iachar(partitionning_strategy(i:i)) - iachar('0')
end do
endif
-
+
case(3)
scotch_strategy = trim(partitionning_strategy)
-
- case default
+
+ case default
print *, 'Invalid partionning method number.'
print *, 'Partionning method', partitionning_method, 'was requested, but is not available.'
stop
end select
-
-
+
+
! read grid parameters
call read_value_double_precision(IIN,IGNORE_JUNK,xmin)
call read_value_double_precision(IIN,IGNORE_JUNK,xmax)
@@ -252,36 +252,36 @@
if ( ngnod == 9 .and. read_external_mesh ) then
print *, 'Number of control nodes must be equal to four when reading from external mesh.'
print *, 'ngnod = 9 is not yet supported.'
- stop
+ stop
endif
-
+
call read_value_logical(IIN,IGNORE_JUNK,initialfield)
call read_value_logical(IIN,IGNORE_JUNK,assign_external_model)
call read_value_logical(IIN,IGNORE_JUNK,TURN_ANISOTROPY_ON)
call read_value_logical(IIN,IGNORE_JUNK,TURN_ATTENUATION_ON)
-
+
if ( read_external_mesh ) then
call read_mesh(mesh_file, nelmnts, elmnts, nnodes, num_start)
-
+
else
! get interface data from external file to count the spectral elements along Z
print *,'Reading interface data from file DATA/',interfacesfile(1:len_trim(interfacesfile)),' to count the spectral elements'
open(unit=IIN_INTERFACES,file='DATA/'//interfacesfile,status='old')
-
+
max_npoints_interface = -1
-
+
! read number of interfaces
call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,number_of_interfaces)
if(number_of_interfaces < 2) stop 'not enough interfaces (minimum is 2)'
! loop on all the interfaces
do interface_current = 1,number_of_interfaces
-
+
call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,npoints_interface_bottom)
if(npoints_interface_bottom < 2) stop 'not enough interface points (minimum is 2)'
max_npoints_interface = max(npoints_interface_bottom,max_npoints_interface)
print *,'Reading ',npoints_interface_bottom,' points for interface ',interface_current
-
+
! loop on all the points describing this interface
xinterface_dummy_previous = -HUGEVAL
do ipoint_current = 1,npoints_interface_bottom
@@ -290,43 +290,43 @@
stop 'interface points must be sorted in increasing X'
xinterface_dummy_previous = xinterface_dummy
enddo
-
+
enddo
-
+
! define number of layers
number_of_layers = number_of_interfaces - 1
-
+
allocate(nz_layer(number_of_layers))
-
+
! loop on all the layers
do ilayer = 1,number_of_layers
-
+
! read number of spectral elements in vertical direction in this layer
call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,nz_layer(ilayer))
if(nz_layer(ilayer) < 1) stop 'not enough spectral elements along Z in layer (minimum is 1)'
print *,'There are ',nz_layer(ilayer),' spectral elements along Z in layer ',ilayer
-
+
enddo
close(IIN_INTERFACES)
-
+
! compute total number of spectral elements in vertical direction
nz = sum(nz_layer)
-
+
print *
print *,'Total number of spectral elements along Z = ',nz
print *
-
+
nxread = nx
nzread = nz
-
+
! multiply by 2 if elements have 9 nodes
if(ngnod == 9) then
nx = nx * 2
nz = nz * 2
nz_layer(:) = nz_layer(:) * 2
endif
-
+
nelmnts = nxread * nzread
allocate(elmnts(0:ngnod*nelmnts-1))
if ( ngnod == 4 ) then
@@ -347,7 +347,7 @@
elmnts(num_elmnt*ngnod) = (j-1)*(nxread+1) + (i-1)
elmnts(num_elmnt*ngnod+1) = (j-1)*(nxread+1) + (i-1) + 1
elmnts(num_elmnt*ngnod+2) = j*(nxread+1) + (i-1) + 1
- elmnts(num_elmnt*ngnod+3) = j*(nxread+1) + (i-1)
+ elmnts(num_elmnt*ngnod+3) = j*(nxread+1) + (i-1)
elmnts(num_elmnt*ngnod+4) = (nxread+1)*(nzread+1) + (j-1)*nxread + (i-1)
elmnts(num_elmnt*ngnod+5) = (nxread+1)*(nzread+1) + nxread*(nzread+1) + (j-1)*(nxread*2+1) + (i-1)*2 + 2
elmnts(num_elmnt*ngnod+6) = (nxread+1)*(nzread+1) + j*nxread + (i-1)
@@ -356,7 +356,7 @@
num_elmnt = num_elmnt + 1
end do
end do
-
+
endif
endif
@@ -413,12 +413,12 @@
print *,'Mxz of the source if moment tensor = ',Mxz
print *,'Multiplying factor = ',factor
-! read constants for attenuation
+! read constants for attenuation
call read_value_integer(IIN,IGNORE_JUNK,N_SLS)
call read_value_double_precision(IIN,IGNORE_JUNK,Qp_attenuation)
call read_value_double_precision(IIN,IGNORE_JUNK,Qs_attenuation)
call read_value_double_precision(IIN,IGNORE_JUNK,f0_attenuation)
-
+
! if source is not a Dirac or Heavyside then f0_attenuation is f0
if(.not. (time_function_type == 4 .or. time_function_type == 5)) then
f0_attenuation = f0
@@ -507,7 +507,7 @@
aniso3(i) = aniso3read
aniso4(i) = aniso4read
enddo
-
+
print *
print *, 'Nb of solid or fluid materials = ',nb_materials
print *
@@ -527,7 +527,7 @@
print *
enddo
-
+
if ( read_external_mesh ) then
call read_mat(materials_file, nelmnts, num_material)
else
@@ -585,7 +585,7 @@
num_material((j-1)*nxread+i) = imaterial_number
enddo
enddo
-
+
enddo
if(minval(num_material) <= 0) stop 'Velocity model not entirely set...'
@@ -711,7 +711,7 @@
enddo
close(IIN_INTERFACES)
-
+
nnodes = (nz+1)*(nx+1)
allocate(nodes_coords(2,nnodes))
if ( ngnod == 4 ) then
@@ -720,34 +720,34 @@
num_node = num_4(i,j,nxread)
nodes_coords(1, num_node) = x(i,j)
nodes_coords(2, num_node) = z(i,j)
-
+
end do
end do
-
+
else
do j = 0, nz
do i = 0, nx
num_node = num_9(i,j,nxread,nzread)
nodes_coords(1, num_node) = x(i,j)
nodes_coords(2, num_node) = z(i,j)
-
+
end do
end do
-
+
endif
else
call read_nodes_coords(nodes_coords_file, nnodes, nodes_coords)
endif
-
+
if ( read_external_mesh ) then
call read_acoustic_surface(free_surface_file, nelem_acoustic_surface, acoustic_surface, &
nelmnts, num_material, ANISOTROPIC_MATERIAL, nb_materials, icodemat, cs, num_start)
-
+
if ( any_abs ) then
call read_abs_surface(absorbing_surface_file, nelemabs, abs_surface, num_start)
endif
-
+
else
! count the number of acoustic free-surface elements
@@ -763,9 +763,9 @@
nelem_acoustic_surface = nelem_acoustic_surface + 1
endif
enddo
-
+
allocate(acoustic_surface(4,nelem_acoustic_surface))
-
+
nelem_acoustic_surface = 0
j = nzread
do i = 1,nxread
@@ -775,12 +775,12 @@
acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
acoustic_surface(2,nelem_acoustic_surface) = 2
acoustic_surface(3,nelem_acoustic_surface) = elmnts(3+ngnod*((j-1)*nxread+i-1))
- acoustic_surface(4,nelem_acoustic_surface) = elmnts(2+ngnod*((j-1)*nxread+i-1))
+ acoustic_surface(4,nelem_acoustic_surface) = elmnts(2+ngnod*((j-1)*nxread+i-1))
endif
end do
endif
-
+
!
!--- definition of absorbing boundaries
!
@@ -789,9 +789,9 @@
if(abstop) nelemabs = nelemabs + nxread
if(absleft) nelemabs = nelemabs + nzread
if(absright) nelemabs = nelemabs + nzread
-
+
allocate(abs_surface(4,nelemabs))
-
+
! generate the list of absorbing elements
if(nelemabs > 0) then
nelemabs = 0
@@ -816,30 +816,30 @@
nelemabs = nelemabs + 1
abs_surface(1,nelemabs) = inumelem-1
abs_surface(2,nelemabs) = 2
- abs_surface(3,nelemabs) = elmnts(3+ngnod*(inumelem-1))
- abs_surface(4,nelemabs) = elmnts(2+ngnod*(inumelem-1))
+ abs_surface(3,nelemabs) = elmnts(3+ngnod*(inumelem-1))
+ abs_surface(4,nelemabs) = elmnts(2+ngnod*(inumelem-1))
endif
if(absleft .and. ix == 1) then
nelemabs = nelemabs + 1
abs_surface(1,nelemabs) = inumelem-1
abs_surface(2,nelemabs) = 2
- abs_surface(3,nelemabs) = elmnts(0+ngnod*(inumelem-1))
- abs_surface(4,nelemabs) = elmnts(3+ngnod*(inumelem-1))
+ abs_surface(3,nelemabs) = elmnts(0+ngnod*(inumelem-1))
+ abs_surface(4,nelemabs) = elmnts(3+ngnod*(inumelem-1))
endif
end do
end do
endif
-
+
endif
-
-
+
+
! compute min and max of X and Z in the grid
print *
print *,'Min and max value of X in the grid = ',minval(nodes_coords(1,:)),maxval(nodes_coords(1,:))
print *,'Min and max value of Z in the grid = ',minval(nodes_coords(2,:)),maxval(nodes_coords(2,:))
print *
-
-
+
+
! ***
! *** create a Gnuplot file that displays the grid
! ***
@@ -902,8 +902,8 @@
print *,'Grid saved in Gnuplot format...'
print *
endif
-
+
!*****************************
! Partitionning
!*****************************
@@ -916,25 +916,25 @@
do i = 0, nelmnts-1
elmnts_bis(i*esize:i*esize+esize-1) = elmnts(i*ngnod:i*ngnod+esize-1)
end do
-
+
if ( nproc > 1 ) then
call mesh2dual_ncommonnodes(nelmnts, (nxread+1)*(nzread+1), elmnts_bis, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
endif
-
+
else
if ( nproc > 1 ) then
call mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
endif
-
+
endif
-
-
+
+
if ( nproc == 1 ) then
part(:) = 0
else
nb_edges = xadj(nelmnts)
-
+
! giving weight to edges and vertices. Currently not used.
call read_weights(nelmnts, vwgt, nb_edges, adjwgt)
@@ -945,7 +945,7 @@
part(iproc*floor(real(nelmnts)/real(nproc)):(iproc+1)*floor(real(nelmnts)/real(nproc))-1) = iproc
end do
part(floor(real(nelmnts)/real(nproc))*(nproc-1):nelmnts-1) = nproc - 1
-
+
case(2)
#ifdef USE_METIS
call Part_metis(nelmnts, xadj, adjncy, vwgt, adjwgt, nproc, nb_edges, edgecut, part, metis_options)
@@ -954,7 +954,7 @@
print *, 'Please recompile with -DUSE_METIS in order to enable use of METIS.'
stop
#endif
-
+
case(3)
#ifdef USE_SCOTCH
call Part_scotch(nelmnts, xadj, adjncy, vwgt, adjwgt, nproc, nb_edges, edgecut, part, scotch_strategy)
@@ -963,11 +963,11 @@
print *, 'Please recompile with -DUSE_SCOTCH in order to enable use of SCOTCH.'
stop
#endif
-
+
end select
-
+
endif
-
+
! beware of fluid solid edges : coupled elements are transfered to the same partition
if ( ngnod == 9 ) then
call acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts_bis, nb_materials, cs, num_material, &
@@ -976,15 +976,15 @@
call acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts, nb_materials, cs, num_material, &
nproc, part, nedges_coupled, edges_coupled)
endif
-
+
! local number of each element for each partition
call Construct_glob2loc_elmnts(nelmnts, part, nproc, glob2loc_elmnts)
-
+
if ( ngnod == 9 ) then
if ( nproc > 1 ) then
deallocate(nnodes_elmnts)
deallocate(nodes_elmnts)
- endif
+ endif
allocate(nnodes_elmnts(0:nnodes-1))
allocate(nodes_elmnts(0:nsize*nnodes-1))
nnodes_elmnts(:) = 0
@@ -992,7 +992,7 @@
do i = 0, ngnod*nelmnts-1
nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/ngnod
nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
-
+
end do
else
if ( nproc < 2 ) then
@@ -1003,14 +1003,14 @@
do i = 0, ngnod*nelmnts-1
nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/ngnod
nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
-
+
end do
- endif
+ endif
endif
-
-
+
+
! local number of each node for each partition
call Construct_glob2loc_nodes(nelmnts, nnodes, nnodes_elmnts, nodes_elmnts, part, nproc, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
@@ -1029,7 +1029,7 @@
allocate(my_nb_interfaces(0:ninterfaces-1))
print *, '05'
endif
-
+
! setting absorbing boundaries by elements instead of edges
if ( any_abs ) then
call merge_abs_boundaries(nelemabs, nelemabs_merge, abs_surface, abs_surface_char, abs_surface_merge, &
@@ -1044,7 +1044,7 @@
! *** generate the databases for the solver
do iproc = 0, nproc-1
-
+
write(prname, "('/Database',i5.5)") iproc
open(unit=15,file='./OUTPUT_FILES'//prname,status='unknown')
@@ -1052,37 +1052,37 @@
write(15,*) '# Database for SPECFEM2D'
write(15,*) '# Dimitri Komatitsch, (c) University of Pau, France'
write(15,*) '#'
-
+
write(15,*) 'Title of the simulation'
write(15,"(a50)") title
-
+
call write_glob2loc_nodes_database(15, iproc, npgeo, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, nnodes, 1)
-
+
call write_partition_database(15, iproc, nspec, nelmnts, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
glob2loc_nodes_parts, glob2loc_nodes, part, num_material, ngnod, 1)
-
+
write(15,*) 'npgeo'
write(15,*) npgeo
-
+
write(15,*) 'gnuplot interpol'
write(15,*) gnuplot,interpol
-
+
write(15,*) 'NTSTEP_BETWEEN_OUTPUT_INFO'
write(15,*) NTSTEP_BETWEEN_OUTPUT_INFO
-
+
write(15,*) 'output_postscript_snapshot output_color_image colors numbers'
write(15,*) output_postscript_snapshot,output_color_image,' 1 0'
write(15,*) 'meshvect modelvect boundvect cutsnaps subsamp sizemax_arrows'
write(15,*) meshvect,modelvect,boundvect,cutsnaps,subsamp,sizemax_arrows
-
+
write(15,*) 'anglerec'
write(15,*) anglerec
-
+
write(15,*) 'initialfield'
write(15,*) initialfield
@@ -1091,27 +1091,27 @@
write(15,*) 'assign_external_model outputgrid TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
write(15,*) assign_external_model,outputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
-
+
write(15,*) 'nt deltat'
write(15,*) nt,deltat
write(15,*) 'source'
write(15,*) source_type,time_function_type,xs,zs,f0,t0,factor,angleforce,Mxx,Mzz,Mxz
-
+
write(15,*) 'attenuation'
write(15,*) N_SLS, Qp_attenuation, Qs_attenuation, f0_attenuation
write(15,*) 'Coordinates of macrobloc mesh (coorg):'
-
+
call write_glob2loc_nodes_database(15, iproc, npgeo, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, nnodes, 2)
-
+
write(15,*) 'numat ngnod nspec pointsdisp plot_lowerleft_corner_only'
write(15,*) nb_materials,ngnod,nspec,pointsdisp,plot_lowerleft_corner_only
-
+
if ( any_abs ) then
call write_abs_merge_database(15, nelemabs_merge, nelemabs_loc, &
abs_surface_char, abs_surface_merge, &
@@ -1121,38 +1121,38 @@
else
nelemabs_loc = 0
endif
-
+
call Write_surface_database(15, nelem_acoustic_surface, acoustic_surface, nelem_acoustic_surface_loc, &
iproc, glob2loc_elmnts, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part, 1)
-
+
call write_fluidsolid_edges_database(15, nedges_coupled, nedges_coupled_loc, &
edges_coupled, glob2loc_elmnts, part, iproc, 1)
-
+
write(15,*) 'nelemabs nelem_acoustic_surface num_fluid_solid_edges'
write(15,*) nelemabs_loc,nelem_acoustic_surface_loc,nedges_coupled_loc
-
-
+
+
write(15,*) 'Material sets (num 1 rho vp vs 0 0) or (num 2 rho c11 c13 c33 c44)'
do i=1,nb_materials
write(15,*) i,icodemat(i),rho(i),cp(i),cs(i),aniso3(i),aniso4(i)
enddo
-
+
write(15,*) 'Arrays kmato and knods for each bloc:'
-
+
call write_partition_database(15, iproc, nspec, nelmnts, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
glob2loc_nodes_parts, glob2loc_nodes, part, num_material, ngnod, 2)
-
- if ( nproc /= 1 ) then
+
+ if ( nproc /= 1 ) then
call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nproc, iproc, ninterfaces, &
my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, 1)
-
+
write(15,*) 'Interfaces:'
write(15,*) my_ninterface, maxval(my_nb_interfaces)
-
+
call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nproc, iproc, ninterfaces, &
my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, 2)
@@ -1160,10 +1160,10 @@
else
write(15,*) 'Interfaces:'
write(15,*) 0, 0
-
+
endif
-
+
write(15,*) 'List of absorbing elements (bottom right top left):'
if ( any_abs ) then
call write_abs_merge_database(15, nelemabs_merge, nelemabs_loc, &
@@ -1172,19 +1172,19 @@
jbegin_left,jend_left,jbegin_right,jend_right, &
glob2loc_elmnts, part, iproc, 2)
endif
-
+
write(15,*) 'List of acoustic free-surface elements:'
call Write_surface_database(15, nelem_acoustic_surface, acoustic_surface, nelem_acoustic_surface_loc, &
iproc, glob2loc_elmnts, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part, 2)
-
+
write(15,*) 'List of acoustic elastic coupled edges:'
call write_fluidsolid_edges_database(15, nedges_coupled, nedges_coupled_loc, &
edges_coupled, glob2loc_elmnts, part, iproc, 2)
end do
-
-
+
+
! print position of the source
print *
print *,'Position (x,z) of the source = ',xs,zs
@@ -1246,8 +1246,8 @@
endif
print *
-
+
end program meshfem2D
! *******************
@@ -1260,45 +1260,45 @@
integer function num(i,j,nx)
implicit none
-
+
integer i,j,nx
-
+
num = j*(nx+1) + i + 1
-
+
end function num
-
+
!--- global node number (when ngnod==4).
integer function num_4(i,j,nx)
implicit none
-
+
integer i,j,nx
-
+
num_4 = j*(nx+1) + i + 1
-
+
end function num_4
-
+
!--- global node number (when ngnod==9).
integer function num_9(i,j,nx,nz)
-
+
implicit none
-
+
integer i,j,nx,nz
-
-
+
+
if ( (mod(i,2) == 0) .and. (mod(j,2) == 0) ) then
num_9 = j/2 * (nx+1) + i/2 + 1
- else
+ else
if ( mod(j,2) == 0 ) then
- num_9 = (nx+1)*(nz+1) + j/2 * nx + ceiling(real(i)/real(2))
- else
+ num_9 = (nx+1)*(nz+1) + j/2 * nx + ceiling(real(i)/real(2))
+ else
num_9 = (nx+1)*(nz+1) + nx*(nz+1) + floor(real(j)/real(2))*(nx*2+1) + i + 1
-
+
endif
endif
-
+
end function num_9
Modified: seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -1,3 +1,4 @@
+
!=====================================================================
!
! Routines from "Numerical Recipes: the Art of Scientific Computing"
Modified: seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -1,13 +1,24 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.2
+! ------------------------------
+!
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
+! University of Pau, France
+!
+! (c) November 2007
+!
+!========================================================================
+
module part_unstruct
implicit none
-
+
include './constants_unstruct.h'
-
contains
-
!-----------------------------------------------
! Read the mesh and storing it in array 'elmnts' (which is allocated here).
! 'num_start' is used to have the numbering of the nodes starting at '0'.
@@ -20,7 +31,7 @@
integer, intent(out) :: nnodes
integer, dimension(:), pointer :: elmnts
integer, intent(out) :: num_start
-
+
integer :: i
print *, trim(filename)
@@ -29,7 +40,7 @@
read(990,*) nelmnts
allocate(elmnts(0:ESIZE*nelmnts-1))
do i = 0, nelmnts-1
- read(990,*) elmnts(i*ESIZE), elmnts(i*ESIZE+1), elmnts(i*ESIZE+2), elmnts(i*ESIZE+3)
+ read(990,*) elmnts(i*ESIZE), elmnts(i*ESIZE+1), elmnts(i*ESIZE+2), elmnts(i*ESIZE+3)
end do
close(990)
@@ -46,7 +57,7 @@
! Read the nodes coordinates and storing it in array 'nodes_coords'
!-----------------------------------------------
subroutine read_nodes_coords(filename, nnodes, nodes_coords)
-
+
character(len=256), intent(in) :: filename
integer, intent(out) :: nnodes
double precision, dimension(:,:), pointer :: nodes_coords
@@ -71,7 +82,7 @@
! Read the material for each element and storing it in array 'num_materials'
!-----------------------------------------------
subroutine read_mat(filename, nelmnts, num_material)
-
+
character(len=256), intent(in) :: filename
integer, intent(in) :: nelmnts
integer, dimension(1:nelmnts), intent(out) :: num_material
@@ -98,12 +109,12 @@
!-----------------------------------------------
subroutine read_acoustic_surface(filename, nelem_acoustic_surface, acoustic_surface, &
nelmnts, num_material, ANISOTROPIC_MATERIAL, nb_materials, icodemat, cs, num_start)
-
+
include './constants.h'
character(len=256), intent(in) :: filename
integer, intent(out) :: nelem_acoustic_surface
- integer, dimension(:,:), pointer :: acoustic_surface
+ integer, dimension(:,:), pointer :: acoustic_surface
integer, intent(in) :: nelmnts
integer, dimension(0:nelmnts-1) :: num_material
integer, intent(in) :: ANISOTROPIC_MATERIAL
@@ -112,39 +123,39 @@
double precision, dimension(1:nb_materials), intent(in) :: cs
integer, intent(in) :: num_start
-
- integer, dimension(:,:), allocatable :: acoustic_surface_tmp
+
+ integer, dimension(:,:), allocatable :: acoustic_surface_tmp
integer :: nelmnts_surface
integer :: i
integer :: imaterial_number
-
-
+
+
open(unit=993, file=trim(filename), form='formatted' , status='old', action='read')
read(993,*) nelmnts_surface
-
+
allocate(acoustic_surface_tmp(4,nelmnts_surface))
-
+
do i = 1, nelmnts_surface
read(993,*) acoustic_surface_tmp(1,i), acoustic_surface_tmp(2,i), acoustic_surface_tmp(3,i), acoustic_surface_tmp(4,i)
-
- end do
+ end do
+
close(993)
acoustic_surface_tmp(1,:) = acoustic_surface_tmp(1,:) - num_start
acoustic_surface_tmp(3,:) = acoustic_surface_tmp(3,:) - num_start
acoustic_surface_tmp(4,:) = acoustic_surface_tmp(4,:) - num_start
-
+
nelem_acoustic_surface = 0
do i = 1, nelmnts_surface
imaterial_number = num_material(acoustic_surface_tmp(1,i))
if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. cs(imaterial_number) < TINYVAL ) then
nelem_acoustic_surface = nelem_acoustic_surface + 1
-
+
end if
end do
-
+
allocate(acoustic_surface(4,nelem_acoustic_surface))
-
+
nelem_acoustic_surface = 0
do i = 1, nelmnts_surface
imaterial_number = num_material(acoustic_surface_tmp(1,i))
@@ -153,41 +164,41 @@
acoustic_surface(:,nelem_acoustic_surface) = acoustic_surface_tmp(:,i)
end if
end do
-
-
+
+
end subroutine read_acoustic_surface
-
-
+
+
!-----------------------------------------------
! Read absorbing surface.
! 'abs_surface' contains 1/ element number, 2/ number of nodes that form the abs surface,
! 3/ first node on the abs surface, 4/ second node on the abs surface, if relevant (if 2/ is equal to 2)
- !-----------------------------------------------
+ !-----------------------------------------------
subroutine read_abs_surface(filename, nelemabs, abs_surface, num_start)
-
+
include './constants.h'
character(len=256), intent(in) :: filename
integer, intent(out) :: nelemabs
- integer, dimension(:,:), pointer :: abs_surface
+ integer, dimension(:,:), pointer :: abs_surface
integer, intent(in) :: num_start
integer :: i
-
-
+
+
open(unit=994, file=trim(filename), form='formatted' , status='old', action='read')
read(994,*) nelemabs
-
+
allocate(abs_surface(4,nelemabs))
-
+
do i = 1, nelemabs
read(994,*) abs_surface(1,i), abs_surface(2,i), abs_surface(3,i), abs_surface(4,i)
-
+
end do
-
+
close(994)
-
+
abs_surface(1,:) = abs_surface(1,:) - num_start
abs_surface(3,:) = abs_surface(3,:) - num_start
abs_surface(4,:) = abs_surface(4,:) - num_start
@@ -196,7 +207,6 @@
end subroutine read_abs_surface
-
!-----------------------------------------------
! Creating dual graph (adjacency is defined by 'ncommonnodes' between two elements).
!-----------------------------------------------
@@ -204,7 +214,7 @@
integer, intent(in) :: nelmnts
integer, intent(in) :: nnodes
- integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
+ integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
integer, dimension(:), pointer :: xadj
integer, dimension(:), pointer :: adjncy
integer, dimension(:), pointer :: nnodes_elmnts
@@ -212,12 +222,12 @@
integer, intent(in) :: ncommonnodes
integer :: i, j, k, l, m, nb_edges
- logical :: is_neighbour
+ logical :: is_neighbour
integer :: num_node, n
integer :: elem_base, elem_target
integer :: connectivity
-
+
allocate(xadj(0:nelmnts))
xadj(:) = 0
allocate(adjncy(0:max_neighbour*nelmnts-1))
@@ -254,16 +264,16 @@
end if
end do
end do
-
+
if ( connectivity >= ncommonnodes) then
-
+
is_neighbour = .false.
-
+
do m = 0, xadj(nodes_elmnts(k+j*nsize))
if ( .not.is_neighbour ) then
if ( adjncy(nodes_elmnts(k+j*nsize)*max_neighbour+m) == nodes_elmnts(l+j*nsize) ) then
is_neighbour = .true.
-
+
end if
end if
end do
@@ -287,31 +297,28 @@
nb_edges = nb_edges + 1
end do
end do
-
+
xadj(nelmnts) = nb_edges
end subroutine mesh2dual_ncommonnodes
-
!-----------------------------------------------
! Read the weight for each vertices and edges of the graph (not curretly used)
!-----------------------------------------------
subroutine read_weights(nelmnts, vwgt, nb_edges, adjwgt)
-
- integer, intent(in) :: nelmnts, nb_edges
+
+ integer, intent(in) :: nelmnts, nb_edges
integer, dimension(:), pointer :: vwgt, adjwgt
-
+
allocate(vwgt(0:nelmnts-1))
allocate(adjwgt(0:nb_edges-1))
-
+
vwgt(:) = 1
adjwgt(:) = 1
-
-
+
end subroutine read_weights
-
!--------------------------------------------------
@@ -345,7 +352,6 @@
end subroutine Construct_glob2loc_elmnts
-
!--------------------------------------------------
! construct local numbering for the nodes in each partition
!--------------------------------------------------
@@ -367,7 +373,6 @@
integer, dimension(0:nparts-1) :: parts_node
integer, dimension(0:nparts-1) :: num_parts
-
allocate(glob2loc_nodes_nparts(0:nnodes))
size_glob2loc_nodes = 0
@@ -392,9 +397,8 @@
end do
+ glob2loc_nodes_nparts(nnodes) = size_glob2loc_nodes
- glob2loc_nodes_nparts(nnodes) = size_glob2loc_nodes
-
allocate(glob2loc_nodes_parts(0:glob2loc_nodes_nparts(nnodes)-1))
allocate(glob2loc_nodes(0:glob2loc_nodes_nparts(nnodes)-1))
@@ -427,14 +431,14 @@
end subroutine Construct_glob2loc_nodes
-
!--------------------------------------------------
! Construct interfaces between each partitions.
- ! Two adjacent elements in distinct partitions make an entry in array tab_interfaces :
- ! 1/ first element, 2/ second element, 3/ number of common nodes, 4/ first node,
+ ! Two adjacent elements in distinct partitions make an entry in array tab_interfaces :
+ ! 1/ first element, 2/ second element, 3/ number of common nodes, 4/ first node,
! 5/ second node, if relevant.
! No interface between acoustic and elastic elements.
!--------------------------------------------------
+
subroutine Construct_interfaces(nelmnts, nparts, part, elmnts, xadj, adjncy, tab_interfaces, &
tab_size_interfaces, ninterfaces, nb_materials, cs_material, num_material)
@@ -450,14 +454,13 @@
integer, dimension(1:nelmnts), intent(in) :: num_material
double precision, dimension(1:nb_materials), intent(in) :: cs_material
integer, intent(in) :: nb_materials
-
+
integer :: num_part, num_part_bis, el, el_adj, num_interface, num_edge, ncommon_nodes, &
num_node, num_node_bis
integer :: i, j
logical :: is_acoustic_el, is_acoustic_el_adj
-
ninterfaces = 0
do i = 0, nparts-1
do j = i+1, nparts-1
@@ -493,7 +496,7 @@
end do
end if
end do
- tab_size_interfaces(num_interface+1) = tab_size_interfaces(num_interface) + num_edge
+ tab_size_interfaces(num_interface+1) = tab_size_interfaces(num_interface) + num_edge
num_edge = 0
num_interface = num_interface + 1
@@ -536,7 +539,7 @@
end do
if ( ncommon_nodes > 0 ) then
tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+2) = ncommon_nodes
- else
+ else
print *, "Error while building interfaces!", ncommon_nodes
end if
num_edge = num_edge + 1
@@ -557,30 +560,30 @@
!--------------------------------------------------
! Write nodes (their coordinates) pertaining to iproc partition in the corresponding Database
!--------------------------------------------------
- subroutine write_glob2loc_nodes_database(IIN_database, iproc, npgeo, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
+ subroutine write_glob2loc_nodes_database(IIN_database, iproc, npgeo, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, nnodes, num_phase)
-
+
integer, intent(in) :: IIN_database
integer, intent(in) :: nnodes, iproc, num_phase
integer, intent(inout) :: npgeo
-
+
double precision, dimension(2,nnodes) :: nodes_coords
integer, dimension(:), pointer :: glob2loc_nodes_nparts
integer, dimension(:), pointer :: glob2loc_nodes_parts
integer, dimension(:), pointer :: glob2loc_nodes
-
+
integer :: i, j
- if ( num_phase == 1 ) then
+ if ( num_phase == 1 ) then
npgeo = 0
-
+
do i = 0, nnodes-1
do j = glob2loc_nodes_nparts(i), glob2loc_nodes_nparts(i+1)-1
if ( glob2loc_nodes_parts(j) == iproc ) then
npgeo = npgeo + 1
-
+
end if
-
+
end do
end do
else
@@ -592,11 +595,10 @@
end do
end do
end if
-
+
end subroutine Write_glob2loc_nodes_database
-
!--------------------------------------------------
! Write elements (their nodes) pertaining to iproc partition in the corresponding Database
!--------------------------------------------------
@@ -615,31 +617,30 @@
integer :: i,j,k
integer, dimension(0:ngnod-1) :: loc_nodes
-
- if ( num_phase == 1 ) then
+ if ( num_phase == 1 ) then
nspec = 0
do i = 0, nelmnts-1
if ( part(i) == iproc ) then
nspec = nspec + 1
-
+
end if
end do
-
- else
+
+ else
do i = 0, nelmnts-1
if ( part(i) == iproc ) then
-
+
do j = 0, ngnod-1
do k = glob2loc_nodes_nparts(elmnts(i*ngnod+j)), glob2loc_nodes_nparts(elmnts(i*ngnod+j)+1)-1
-
+
if ( glob2loc_nodes_parts(k) == iproc ) then
loc_nodes(j) = glob2loc_nodes(k)
-
+
end if
end do
-
+
end do
write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(i+1), (loc_nodes(k)+1, k=0,ngnod-1)
end if
@@ -650,15 +651,13 @@
end subroutine write_partition_database
-
!--------------------------------------------------
! Write interfaces (element and common nodes) pertaining to iproc partition in the corresponding Database
!--------------------------------------------------
subroutine Write_interfaces_database(IIN_database, tab_interfaces, tab_size_interfaces, nparts, iproc, ninterfaces, &
my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, num_phase)
-
-
+
integer, intent(in) :: IIN_database
integer, intent(in) :: iproc
integer, intent(in) :: nparts
@@ -672,23 +671,21 @@
integer, dimension(:), pointer :: glob2loc_nodes_nparts
integer, dimension(:), pointer :: glob2loc_nodes_parts
integer, dimension(:), pointer :: glob2loc_nodes
-
integer, dimension(2) :: local_nodes
integer :: local_elmnt
integer :: num_phase
-
+
integer :: i, j, k, l
integer :: num_interface
-
num_interface = 0
-
- if ( num_phase == 1 ) then
-
+
+ if ( num_phase == 1 ) then
+
my_interfaces(:) = 0
my_nb_interfaces(:) = 0
-
+
do i = 0, nparts-1
do j = i+1, nparts-1
if ( (tab_size_interfaces(num_interface) < tab_size_interfaces(num_interface+1)) .and. &
@@ -700,33 +697,33 @@
end do
end do
my_ninterface = sum(my_interfaces(:))
-
+
else
-
+
do i = 0, nparts-1
do j = i+1, nparts-1
if ( my_interfaces(num_interface) == 1 ) then
if ( i == iproc ) then
write(IIN_database,*) j, my_nb_interfaces(num_interface)
- else
+ else
write(IIN_database,*) i, my_nb_interfaces(num_interface)
end if
-
+
do k = tab_size_interfaces(num_interface), tab_size_interfaces(num_interface+1)-1
if ( i == iproc ) then
local_elmnt = glob2loc_elmnts(tab_interfaces(k*5+0))+1
- else
+ else
local_elmnt = glob2loc_elmnts(tab_interfaces(k*5+1))+1
end if
- if ( tab_interfaces(k*5+2) == 1 ) then
+ if ( tab_interfaces(k*5+2) == 1 ) then
do l = glob2loc_nodes_nparts(tab_interfaces(k*5+3)), &
glob2loc_nodes_nparts(tab_interfaces(k*5+3)+1)-1
if ( glob2loc_nodes_parts(l) == iproc ) then
local_nodes(1) = glob2loc_nodes(l)+1
end if
end do
-
+
write(IIN_database,*) local_elmnt, tab_interfaces(k*5+2), local_nodes(1), -1
else
if ( tab_interfaces(k*5+2) == 2 ) then
@@ -743,75 +740,70 @@
end if
end do
write(IIN_database,*) local_elmnt, tab_interfaces(k*5+2), local_nodes(1), local_nodes(2)
- else
+ else
write(IIN_database,*) "erreur_write_interface_", tab_interfaces(k*5+2)
end if
end if
end do
-
+
end if
-
+
num_interface = num_interface + 1
end do
end do
-
+
end if
-
-
+
end subroutine Write_interfaces_database
-
!--------------------------------------------------
! Write a surface (elements and nodes on the surface) pertaining to iproc partition in the corresponding Database
!--------------------------------------------------
+
subroutine Write_surface_database(IIN_database, nsurface, surface, &
nsurface_loc, iproc, glob2loc_elmnts, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part, num_phase)
-
-
+
integer, intent(in) :: IIN_database
integer, intent(in) :: iproc
integer :: nsurface
integer :: nsurface_loc
integer, dimension(:,:), pointer :: surface
-
+
integer, dimension(:), pointer :: glob2loc_elmnts
integer, dimension(:), pointer :: glob2loc_nodes_nparts
integer, dimension(:), pointer :: glob2loc_nodes_parts
integer, dimension(:), pointer :: glob2loc_nodes
integer, dimension(:), pointer :: part
-
-
+
integer, dimension(2) :: local_nodes
integer :: local_elmnt
integer :: num_phase
-
+
integer :: i, l
-
-
- if ( num_phase == 1 ) then
-
+ if ( num_phase == 1 ) then
+
nsurface_loc = 0
-
+
do i = 1, nsurface
if ( part(surface(1,i)) == iproc ) then
nsurface_loc = nsurface_loc + 1
-
+
end if
end do
-
+
else
-
+
nsurface_loc = 0
-
+
do i = 1, nsurface
if ( part(surface(1,i)) == iproc ) then
nsurface_loc = nsurface_loc + 1
-
+
local_elmnt = glob2loc_elmnts(surface(1,i)) + 1
-
+
if ( surface(2,i) == 1 ) then
do l = glob2loc_nodes_nparts(surface(3,i)), &
glob2loc_nodes_nparts(surface(3,i)+1)-1
@@ -819,7 +811,7 @@
local_nodes(1) = glob2loc_nodes(l)+1
end if
end do
-
+
write(IIN_database,*) local_elmnt, surface(2,i), local_nodes(1), -1
end if
if ( surface(2,i) == 2 ) then
@@ -835,27 +827,26 @@
local_nodes(2) = glob2loc_nodes(l)+1
end if
end do
-
+
write(IIN_database,*) local_elmnt, surface(2,i), local_nodes(1), local_nodes(2)
end if
-
+
end if
-
+
end do
-
+
end if
-
-
+
end subroutine Write_surface_database
-
!--------------------------------------------------
! Set absorbing boundaries by elements instead of edges.
- ! Excludes points that have both absorbing condition and coupled fluid/solid relation (this is the
+ ! Excludes points that have both absorbing condition and coupled fluid/solid relation (this is the
! reason arrays ibegin_..., iend_... were included here).
! Under development : exluding points that have two different normal.
!--------------------------------------------------
+
subroutine merge_abs_boundaries(nelemabs, nelemabs_merge, abs_surface, abs_surface_char, abs_surface_merge, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right, &
@@ -866,7 +857,6 @@
implicit none
include 'constants.h'
-
integer, intent(inout) :: nelemabs
integer, intent(out) :: nelemabs_merge
integer, dimension(:,:), pointer :: abs_surface
@@ -891,19 +881,18 @@
integer :: i
integer :: temp
integer :: iedge, inode1, inode2
-
-
-
+
+
allocate(abs_surface_char(4,nelemabs))
allocate(abs_surface_merge(nelemabs))
abs_surface_char(:,:) = .false.
abs_surface_merge(:) = -1
-
+
nedge_bound = nelemabs
nb_elmnts_abs = 0
do num_edge = 1, nedge_bound
-
+
match = 0
do i = 1, nb_elmnts_abs
if ( abs_surface(1,num_edge) == abs_surface_merge(i) ) then
@@ -911,82 +900,81 @@
exit
end if
end do
-
+
if ( match == 0 ) then
nb_elmnts_abs = nb_elmnts_abs + 1
match = nb_elmnts_abs
end if
abs_surface_merge(match) = abs_surface(1,num_edge)
-
-
+
+
if ( (abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+0) .and. &
abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+1)) ) then
- abs_surface_char(1,match) = .true.
-
+ abs_surface_char(1,match) = .true.
+
end if
-
+
if ( (abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+0) .and. &
abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+1)) ) then
temp = abs_surface(4,num_edge)
abs_surface(4,num_edge) = abs_surface(3,num_edge)
abs_surface(3,num_edge) = temp
abs_surface_char(1,match) = .true.
-
+
end if
-
+
if ( (abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+0) .and. &
abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+3)) ) then
abs_surface_char(4,match) = .true.
-
+
end if
-
+
if ( (abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+0) .and. &
abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+3)) ) then
temp = abs_surface(4,num_edge)
abs_surface(4,num_edge) = abs_surface(3,num_edge)
abs_surface(3,num_edge) = temp
abs_surface_char(4,match) = .true.
-
+
end if
-
+
if ( (abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+1) .and. &
abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+2)) ) then
abs_surface_char(2,match) = .true.
end if
-
+
if ( (abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+1) .and. &
abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+2)) ) then
temp = abs_surface(4,num_edge)
abs_surface(4,num_edge) = abs_surface(3,num_edge)
abs_surface(3,num_edge) = temp
abs_surface_char(2,match) = .true.
-
+
end if
-
+
if ( (abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+2) .and. &
abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+3)) ) then
temp = abs_surface(4,num_edge)
abs_surface(4,num_edge) = abs_surface(3,num_edge)
abs_surface(3,num_edge) = temp
abs_surface_char(3,match) = .true.
-
+
end if
-
+
if ( (abs_surface(4,num_edge) == elmnts(ngnod*abs_surface_merge(match)+2) .and. &
abs_surface(3,num_edge) == elmnts(ngnod*abs_surface_merge(match)+3)) ) then
abs_surface_char(3,match) = .true.
-
+
end if
-
+
end do
-
+
nelemabs_merge = nb_elmnts_abs
-
-
+
allocate(ibegin_bottom(nelemabs_merge))
- allocate(iend_bottom(nelemabs_merge))
+ allocate(iend_bottom(nelemabs_merge))
allocate(jbegin_right(nelemabs_merge))
allocate(jend_right(nelemabs_merge))
allocate(ibegin_top(nelemabs_merge))
@@ -1008,10 +996,9 @@
if (cs_material(i) < TINYVAL) then
is_acoustic(i) = .true.
end if
-
+
end do
-
do num_edge = 1, nedge_bound
match = 0
@@ -1023,91 +1010,91 @@
end do
if ( is_acoustic(num_material(abs_surface(1,num_edge)+1)) ) then
-
+
do iedge = 1, nedges_coupled
-
+
do inode1 = 0, 3
if ( abs_surface(3,num_edge) == elmnts(ngnod*edges_coupled(1,iedge)+inode1) ) then
do inode2 = 0, 3
if ( abs_surface(3,num_edge) == elmnts(ngnod*edges_coupled(2,iedge)+inode2) ) then
if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+0) .and. &
- abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+1) ) then
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+1) ) then
ibegin_bottom(match) = 2
-
+
end if
if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+1) .and. &
- abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
jbegin_right(match) = 2
-
+
end if
if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+3) .and. &
- abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
ibegin_top(match) = 2
-
+
end if
if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+0) .and. &
- abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+3) ) then
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+3) ) then
jbegin_left(match) = 2
-
+
end if
-
+
end if
end do
-
+
end if
-
+
if ( abs_surface(4,num_edge) == elmnts(ngnod*edges_coupled(1,iedge)+inode1) ) then
do inode2 = 0, 3
if ( abs_surface(4,num_edge) == elmnts(ngnod*edges_coupled(2,iedge)+inode2) ) then
if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+0) .and. &
- abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+1) ) then
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+1) ) then
iend_bottom(match) = NGLLX - 1
-
+
end if
if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+1) .and. &
- abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
jend_right(match) = NGLLZ - 1
-
+
end if
if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+3) .and. &
- abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+2) ) then
iend_top(match) = NGLLX - 1
-
+
end if
if ( abs_surface(3,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+0) .and. &
- abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+3) ) then
+ abs_surface(4,num_edge) == elmnts(ngnod*abs_surface(1,num_edge)+3) ) then
jend_left(match) = NGLLZ - 1
-
+
end if
end if
end do
-
+
end if
-
+
end do
-
-
+
+
end do
-
+
end if
end do
-
end subroutine merge_abs_boundaries
-
-
+
+
!--------------------------------------------------
! Write abs surface (elements and nodes on the surface) pertaining to iproc partition in the corresponding Database
- !--------------------------------------------------
+ !--------------------------------------------------
+
subroutine write_abs_merge_database(IIN_database, nelemabs_merge, nelemabs_loc, &
abs_surface_char, abs_surface_merge, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right, &
glob2loc_elmnts, part, iproc, num_phase)
-
+
implicit none
-
+
integer, intent(in) :: IIN_database
integer, intent(in) :: nelemabs_merge
integer, intent(inout) :: nelemabs_loc
@@ -1119,10 +1106,9 @@
integer, intent(in) :: num_phase
integer, dimension(:), pointer :: ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right
-
+
integer :: i
-
-
+
if ( num_phase == 1 ) then
nelemabs_loc = 0
do i = 1, nelemabs_merge
@@ -1133,31 +1119,30 @@
else
do i = 1, nelemabs_merge
if ( part(abs_surface_merge(i)) == iproc ) then
-
+
write(IIN_database,*) glob2loc_elmnts(abs_surface_merge(i))+1, abs_surface_char(1,i), &
abs_surface_char(2,i), abs_surface_char(3,i), abs_surface_char(4,i), &
- ibegin_bottom(i), iend_bottom(i), &
- jbegin_right(i), jend_right(i), &
- ibegin_top(i), iend_top(i), &
+ ibegin_bottom(i), iend_bottom(i), &
+ jbegin_right(i), jend_right(i), &
+ ibegin_top(i), iend_top(i), &
jbegin_left(i), jend_left(i)
-
+
end if
-
+
end do
end if
-
-
+
+
end subroutine write_abs_merge_database
-
-
-
+
+
#ifdef USE_METIS
!--------------------------------------------------
! Partitioning using METIS
!--------------------------------------------------
subroutine Part_metis(nelmnts, xadj, adjncy, vwgt, adjwgt, nparts, nb_edges, edgecut, part, metis_options)
- integer, intent(in) :: nelmnts, nparts, nb_edges
+ integer, intent(in) :: nelmnts, nparts, nb_edges
integer, intent(inout) :: edgecut
integer, dimension(0:nelmnts), intent(in) :: xadj
integer, dimension(0:max_neighbour*nelmnts-1), intent(in) :: adjncy
@@ -1184,16 +1169,15 @@
#endif
-
#ifdef USE_SCOTCH
!--------------------------------------------------
! Partitioning using SCOTCH
!--------------------------------------------------
subroutine Part_scotch(nelmnts, xadj, adjncy, vwgt, adjwgt, nparts, nedges, edgecut, part, scotch_strategy)
-
+
include 'scotchf.h'
- integer, intent(in) :: nelmnts, nparts, nedges
+ integer, intent(in) :: nelmnts, nparts, nedges
integer, intent(inout) :: edgecut
integer, dimension(0:nelmnts), intent(in) :: xadj
integer, dimension(0:max_neighbour*nelmnts-1), intent(in) :: adjncy
@@ -1206,7 +1190,6 @@
character(len=256), intent(in) :: scotch_strategy
integer :: IERR
-
edgecut = vwgt(0)
edgecut = 0
@@ -1215,13 +1198,13 @@
PRINT *, 'ERROR : MAIN : Cannot initialize strat'
STOP
ENDIF
-
+
call scotchfstratgraphmap (SCOTCHSTRAT(1), trim(scotch_strategy), IERR)
IF (IERR .NE. 0) THEN
PRINT *, 'ERROR : MAIN : Cannot build strat'
STOP
ENDIF
-
+
CALL SCOTCHFGRAPHINIT (SCOTCHGRAPH (1), IERR)
IF (IERR .NE. 0) THEN
PRINT *, 'ERROR : MAIN : Cannot initialize graph'
@@ -1232,18 +1215,18 @@
xadj (0), xadj (0), &
xadj (0), xadj (0), &
nedges, &
- adjncy (0), adjwgt (0), IERR)
+ adjncy (0), adjwgt (0), IERR)
IF (IERR .NE. 0) THEN
PRINT *, 'ERROR : MAIN : Cannot build graph'
STOP
ENDIF
-
+
CALL SCOTCHFGRAPHCHECK (SCOTCHGRAPH (1), IERR)
IF (IERR .NE. 0) THEN
PRINT *, 'ERROR : MAIN : Invalid check'
STOP
ENDIF
-
+
call scotchfgraphpart (SCOTCHGRAPH (1), nparts, SCOTCHSTRAT(1), part(0), IERR)
IF (IERR .NE. 0) THEN
PRINT *, 'ERROR : MAIN : Cannot part graph'
@@ -1262,19 +1245,20 @@
STOP
ENDIF
-
+
end subroutine Part_scotch
#endif
-
!--------------------------------------------------
! Repartitioning : two coupled acoustic/elastic elements are transfered to the same partition
!--------------------------------------------------
+
subroutine acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts, nb_materials, cs_material, num_material, &
nproc, part, nedges_coupled, edges_coupled)
-
+
implicit none
+
include 'constants.h'
integer, intent(in) :: nelmnts, nnodes, nproc, nb_materials
@@ -1284,29 +1268,26 @@
integer, dimension(:), pointer :: part
integer, intent(out) :: nedges_coupled
integer, dimension(:,:), pointer :: edges_coupled
-
-
+
+
logical, dimension(nb_materials) :: is_acoustic
integer, dimension(:), pointer :: xadj
integer, dimension(:), pointer :: adjncy
integer, dimension(:), pointer :: nodes_elmnts
integer, dimension(:), pointer :: nnodes_elmnts
-
+
integer :: i, num_edge
integer :: el, el_adj
logical :: is_repartitioned
-
-
is_acoustic(:) = .false.
do i = 1, nb_materials
if (cs_material(i) < TINYVAL) then
is_acoustic(i) = .true.
end if
-
+
end do
-
call mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,2)
nedges_coupled = 0
@@ -1316,15 +1297,15 @@
if ( .not. is_acoustic(num_material(adjncy(el_adj)+1)) ) then
nedges_coupled = nedges_coupled + 1
end if
-
+
end do
end if
end do
-
+
print *, 'nedges_coupled', nedges_coupled
allocate(edges_coupled(2,nedges_coupled))
-
+
nedges_coupled = 0
do el = 0, nelmnts-1
if ( is_acoustic(num_material(el+1)) ) then
@@ -1334,11 +1315,11 @@
edges_coupled(1,nedges_coupled) = el
edges_coupled(2,nedges_coupled) = adjncy(el_adj)
end if
-
+
end do
end if
end do
-
+
do i = 1, nedges_coupled*nproc
is_repartitioned = .false.
do num_edge = 1, nedges_coupled
@@ -1350,27 +1331,26 @@
end if
is_repartitioned = .true.
end if
-
+
end do
if ( .not. is_repartitioned ) then
exit
end if
end do
-
end subroutine acoustic_elastic_repartitioning
-
!--------------------------------------------------
- ! Write fluid/solid edges (fluid elements and corresponding solid elements)
+ ! Write fluid/solid edges (fluid elements and corresponding solid elements)
! pertaining to iproc partition in the corresponding Database
- !--------------------------------------------------
+ !--------------------------------------------------
+
subroutine write_fluidsolid_edges_database(IIN_database, nedges_coupled, nedges_coupled_loc, &
edges_coupled, glob2loc_elmnts, part, iproc, num_phase)
-
+
implicit none
-
+
integer, intent(in) :: IIN_database
integer, intent(in) :: nedges_coupled
integer, intent(inout) :: nedges_coupled_loc
@@ -1379,11 +1359,9 @@
integer, dimension(:), pointer :: part
integer, intent(in) :: iproc
integer, intent(in) :: num_phase
-
-
+
integer :: i
-
-
+
if ( num_phase == 1 ) then
nedges_coupled_loc = 0
do i = 1, nedges_coupled
@@ -1394,16 +1372,15 @@
else
do i = 1, nedges_coupled
if ( part(edges_coupled(1,i)) == iproc ) then
- write(IIN_database,*) glob2loc_elmnts(edges_coupled(1,i))+1, glob2loc_elmnts(edges_coupled(2,i))+1
-
+ write(IIN_database,*) glob2loc_elmnts(edges_coupled(1,i))+1, glob2loc_elmnts(edges_coupled(2,i))+1
+
end if
-
+
end do
end if
-
-
-end subroutine write_fluidsolid_edges_database
+end subroutine write_fluidsolid_edges_database
end module part_unstruct
+
Modified: seismo/2D/SPECFEM2D/trunk/plotgll.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/plotgll.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/plotgll.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/plotpost.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/plotpost.F90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/plotpost.F90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -107,7 +107,7 @@
double precision, dimension(:,:), allocatable :: RGB_recv
integer :: nspec_recv
integer :: buffer_offset, RGB_offset
-
+
integer :: nb_coorg_per_elem, nb_color_per_elem
integer :: iproc, num_spec
integer :: ier
@@ -1349,7 +1349,7 @@
write(IOUT,*) 'X min, max = ',xmin,xmax
write(IOUT,*) 'Z min, max = ',zmin,zmax
end if
-
+
! ratio of physical page size/size of the domain meshed
ratio_page = min(rpercentz*sizez/(zmax-zmin),rpercentx*sizex/(xmax-xmin)) / 100.d0
@@ -1593,7 +1593,7 @@
else
buffer_offset = buffer_offset + 1
coorg_send(1,buffer_offset) = xw
- coorg_send(2,buffer_offset) = zw
+ coorg_send(2,buffer_offset) = zw
end if
xw = coord(1,ibool(i+subsamp,j+subsamp,ispec))
@@ -1607,7 +1607,7 @@
else
buffer_offset = buffer_offset + 1
coorg_send(1,buffer_offset) = xw
- coorg_send(2,buffer_offset) = zw
+ coorg_send(2,buffer_offset) = zw
end if
xw = coord(1,ibool(i,j+subsamp,ispec))
@@ -1621,7 +1621,7 @@
else
buffer_offset = buffer_offset + 1
coorg_send(1,buffer_offset) = xw
- coorg_send(2,buffer_offset) = zw
+ coorg_send(2,buffer_offset) = zw
end if
! display P-velocity model using gray levels
@@ -1636,10 +1636,10 @@
enddo
enddo
-
+
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
allocate(coorg_recv(2,nspec_recv*((NGLLX-subsamp)/subsamp)*((NGLLX-subsamp)/subsamp)*4))
@@ -1648,7 +1648,7 @@
MPI_DOUBLE_PRECISION, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
call MPI_RECV (RGB_recv(1,1), nspec_recv*((NGLLX-subsamp)/subsamp)*((NGLLX-subsamp)/subsamp), &
MPI_DOUBLE_PRECISION, iproc, 42, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
buffer_offset = 0
RGB_offset = 0
do ispec = 1, nspec_recv
@@ -1678,13 +1678,13 @@
MPI_DOUBLE_PRECISION, 0, 42, MPI_COMM_WORLD, ier)
call MPI_SEND (RGB_send(1,1), nspec*((NGLLX-subsamp)/subsamp)*((NGLLX-subsamp)/subsamp), &
MPI_DOUBLE_PRECISION, 0, 42, MPI_COMM_WORLD, ier)
-
+
deallocate(coorg_send)
deallocate(RGB_send)
-
+
end if
-
+
#endif
@@ -1699,9 +1699,9 @@
write(24,*) '% spectral element mesh'
write(24,*) '%'
end if
-
+
if ( myrank /= 0 ) then
-
+
if ( ngnod == 4 ) then
if ( numbers == 1 ) then
allocate(coorg_send(2,nspec*5))
@@ -1731,7 +1731,7 @@
end if
end if
end if
-
+
end if
buffer_offset = 0
RGB_offset = 0
@@ -1912,7 +1912,7 @@
RGB_offset = RGB_offset + 1
color_send(RGB_offset) = icol
end if
-
+
endif
if ( myrank == 0 ) then
@@ -1934,7 +1934,7 @@
zw = (zw-zmin)*ratio_page + orig_z
xw = xw * centim
zw = zw * centim
-
+
if ( myrank == 0 ) then
if(colors == 1) write(24,*) '1 setgray'
end if
@@ -1959,10 +1959,10 @@
enddo
-
+
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 43, MPI_COMM_WORLD, request_mpi_status, ier)
nb_coorg_per_elem = 1
@@ -1990,7 +1990,7 @@
MPI_DOUBLE_PRECISION, iproc, 43, MPI_COMM_WORLD, request_mpi_status, ier)
call MPI_RECV (color_recv(1), nspec_recv*nb_coorg_per_elem, &
MPI_INTEGER, iproc, 43, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
buffer_offset = 0
RGB_offset = 0
num_spec = nspec
@@ -2027,9 +2027,9 @@
buffer_offset = buffer_offset + 1
write(24,681) coorg_recv(1,buffer_offset), coorg_recv(2,buffer_offset)
end do
-
+
end if
-
+
write(24,*) 'CO'
if ( colors == 1 ) then
if(meshvect) then
@@ -2085,13 +2085,13 @@
call MPI_SEND (color_send(1), nspec*nb_color_per_elem, &
MPI_INTEGER, 0, 43, MPI_COMM_WORLD, ier)
end if
-
+
deallocate(coorg_send)
deallocate(color_send)
end if
-
+
#endif
@@ -2169,18 +2169,18 @@
enddo
end if
-
+
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 44, MPI_COMM_WORLD, request_mpi_status, ier)
if ( nspec_recv > 0 ) then
allocate(coorg_recv(4,nspec_recv))
call MPI_RECV (coorg_recv(1,1), 4*nspec_recv, &
MPI_DOUBLE_PRECISION, iproc, 44, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
buffer_offset = 0
do ispec = 1, nspec_recv
buffer_offset = buffer_offset + 1
@@ -2197,12 +2197,12 @@
MPI_DOUBLE_PRECISION, 0, 44, MPI_COMM_WORLD, ier)
deallocate(coorg_send)
end if
-
+
end if
-
-#endif
+#endif
+
if ( myrank == 0 ) then
write(24,*) '0 setgray'
write(24,*) '0.01 CM setlinewidth'
@@ -2257,18 +2257,18 @@
enddo
end if
-
+
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 44, MPI_COMM_WORLD, request_mpi_status, ier)
if ( nspec_recv > 0 ) then
allocate(coorg_recv(4,nspec_recv))
call MPI_RECV (coorg_recv(1,1), 4*nspec_recv, &
MPI_DOUBLE_PRECISION, iproc, 44, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
buffer_offset = 0
do ispec = 1, nspec_recv
buffer_offset = buffer_offset + 1
@@ -2285,19 +2285,19 @@
MPI_DOUBLE_PRECISION, 0, 44, MPI_COMM_WORLD, ier)
deallocate(coorg_send)
end if
-
+
end if
-
-#endif
+#endif
+
if ( myrank == 0 ) then
write(24,*) '0 setgray'
write(24,*) '0.01 CM setlinewidth'
end if
-
+
!
!---- draw the fluid-solid coupling edges with a thick color line
!
@@ -2371,17 +2371,17 @@
enddo
-
+
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 45, MPI_COMM_WORLD, request_mpi_status, ier)
if ( nspec_recv > 0 ) then
allocate(coorg_recv(4,nspec_recv))
call MPI_RECV (coorg_recv(1,1), 4*nspec_recv, &
MPI_DOUBLE_PRECISION, iproc, 45, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
buffer_offset = 0
do ispec = 1, nspec_recv
buffer_offset = buffer_offset + 1
@@ -2400,10 +2400,10 @@
deallocate(coorg_send)
end if
end if
-
-#endif
+#endif
+
if ( myrank == 0 ) then
write(24,*) '0 setgray'
write(24,*) '0.01 CM setlinewidth'
@@ -2440,7 +2440,7 @@
if ( myrank == 0 ) then
write(IOUT,*) 'Interpolating the vector field...'
end if
-
+
! option to plot only lowerleft corner value to avoid very large files if dense meshes
if(plot_lowerleft_corner_only) then
pointsdisp_loop = 1
@@ -2450,7 +2450,7 @@
if ( myrank /= 0 ) then
allocate(coorg_send(8,nspec*pointsdisp_loop*pointsdisp_loop))
-
+
end if
buffer_offset = 0
@@ -2540,7 +2540,7 @@
enddo
ch2(index_char) = ch1(line_length)
write(24,"(100(a1))") (ch2(ii), ii=1,index_char)
-
+
else
buffer_offset = buffer_offset + 1
coorg_send(1,buffer_offset) = xb
@@ -2552,7 +2552,7 @@
coorg_send(7,buffer_offset) = x1
coorg_send(8,buffer_offset) = z1
end if
-
+
endif
enddo
@@ -2562,14 +2562,14 @@
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 46, MPI_COMM_WORLD, request_mpi_status, ier)
if ( nspec_recv > 0 ) then
allocate(coorg_recv(8,nspec_recv))
call MPI_RECV (coorg_recv(1,1), 8*nspec_recv, &
MPI_DOUBLE_PRECISION, iproc, 46, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
buffer_offset = 0
do ispec = 1, nspec_recv
buffer_offset = buffer_offset + 1
@@ -2581,7 +2581,7 @@
! suppress leading white spaces again, if any
postscript_line = adjustl(postscript_line)
-
+
line_length = len_trim(postscript_line)
index_char = 1
first = .false.
@@ -2607,18 +2607,18 @@
MPI_DOUBLE_PRECISION, 0, 46, MPI_COMM_WORLD, ier)
deallocate(coorg_send)
end if
-
+
end if
-
-#endif
+#endif
+
! draw the vectors at the nodes of the mesh if we do not interpolate the display on a regular grid
else
if ( myrank /= 0 ) then
allocate(coorg_send(8,npoin))
-
+
end if
buffer_offset = 0
@@ -2701,14 +2701,14 @@
#ifdef USE_MPI
if (myrank == 0 ) then
-
+
do iproc = 1, nproc-1
call MPI_RECV (nspec_recv, 1, MPI_INTEGER, iproc, 47, MPI_COMM_WORLD, request_mpi_status, ier)
if ( nspec_recv > 0 ) then
allocate(coorg_recv(8,nspec_recv))
call MPI_RECV (coorg_recv(1,1), 8*nspec_recv, &
MPI_DOUBLE_PRECISION, iproc, 47, MPI_COMM_WORLD, request_mpi_status, ier)
-
+
buffer_offset = 0
do ispec = 1, nspec_recv
buffer_offset = buffer_offset + 1
@@ -2720,7 +2720,7 @@
! suppress leading white spaces again, if any
postscript_line = adjustl(postscript_line)
-
+
line_length = len_trim(postscript_line)
index_char = 1
first = .false.
@@ -2747,10 +2747,10 @@
deallocate(coorg_send)
end if
end if
-
-#endif
+#endif
+
endif
if ( myrank == 0 ) then
Modified: seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/recompute_jacobian.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/recompute_jacobian.f90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/recompute_jacobian.f90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -379,7 +379,7 @@
read(IIN,"(a80)") datlin
read(IIN,*) NTSTEP_BETWEEN_OUTPUT_INFO
-
+
read(IIN,"(a80)") datlin
read(IIN,*) output_postscript_snapshot,output_color_image,colors,numbers
@@ -416,7 +416,7 @@
endif
NTSTEP_BETWEEN_OUTPUT_SEISMO = min(NSTEP,NTSTEP_BETWEEN_OUTPUT_INFO)
-
+
!
!---- read source information
!
@@ -434,7 +434,7 @@
!
if(.not. initialfield) then
if (source_type == 1) then
- if ( myrank == 0 ) then
+ if ( myrank == 0 ) then
write(IOUT,212) x_source,z_source,f0,t0,factor,angleforce
endif
else if(source_type == 2) then
@@ -1066,7 +1066,7 @@
#ifdef USE_MPI
if ( nproc > 1 ) then
! preparing for MPI communications
- allocate(mask_ispec_inner_outer(nspec))
+ allocate(mask_ispec_inner_outer(nspec))
mask_ispec_inner_outer(:) = .false.
call prepare_assemble_MPI (nspec,ibool, &
@@ -1081,8 +1081,8 @@
mask_ispec_inner_outer &
)
- nspec_outer = count(mask_ispec_inner_outer)
- nspec_inner = nspec - nspec_outer
+ nspec_outer = count(mask_ispec_inner_outer)
+ nspec_inner = nspec - nspec_outer
allocate(ispec_outer_to_glob(nspec_outer))
allocate(ispec_inner_to_glob(nspec_inner))
@@ -1097,8 +1097,8 @@
else
num_ispec_inner = num_ispec_inner + 1
ispec_inner_to_glob(num_ispec_inner) = ispec
-
- endif
+
+ endif
enddo
max_ibool_interfaces_size_ac = maxval(nibool_interfaces_acoustic(:))
@@ -1139,17 +1139,17 @@
ninterface, max_interface_size, max_ibool_interfaces_size_ac, max_ibool_interfaces_size_el, &
ibool_interfaces_acoustic,ibool_interfaces_elastic, nibool_interfaces_acoustic,nibool_interfaces_elastic, my_neighbours)
- else
+ else
ninterface_acoustic = 0
ninterface_elastic = 0
-
+
num_ispec_outer = 0
num_ispec_inner = 0
allocate(mask_ispec_inner_outer(1))
nspec_outer = 0
nspec_inner = nspec
-
+
allocate(ispec_inner_to_glob(nspec_inner))
do ispec = 1, nspec
ispec_inner_to_glob(ispec) = ispec
@@ -1164,7 +1164,7 @@
nspec_outer = 0
nspec_inner = nspec
-
+
allocate(ispec_inner_to_glob(nspec_inner))
do ispec = 1, nspec
ispec_inner_to_glob(ispec) = ispec
@@ -1484,9 +1484,9 @@
close(55)
endif
-! nb_proc_source is the number of processes that own the source (the nearest point). It can be greater
+! nb_proc_source is the number of processes that own the source (the nearest point). It can be greater
! than one if the nearest point is on the interface between several partitions with an explosive source.
-! since source contribution is linear, the source_time_function is cut down by that number (it would have been similar
+! since source contribution is linear, the source_time_function is cut down by that number (it would have been similar
! if we just had elected one of those processes).
source_time_function(:) = source_time_function(:) / nb_proc_source
@@ -1626,7 +1626,7 @@
enddo
enddo
-
+
if ( myrank == 0 ) then
print *,'End of fluid/solid edge detection'
print *
@@ -1734,7 +1734,7 @@
! getting velocity for each local pixels
image_color_vp_display(:,:) = 0.d0
-
+
do k = 1, nb_pixel_loc
j = ceiling(real(num_pixel_loc(k)) / real(NX_IMAGE_color))
i = num_pixel_loc(k) - (j-1)*NX_IMAGE_color
@@ -1776,8 +1776,8 @@
! initialize variables for writing seismograms
seismo_offset = 0
seismo_current = 0
-
+
! *********************************************************
! ************* MAIN LOOP OVER THE TIME STEPS *************
! *********************************************************
@@ -1918,7 +1918,7 @@
hprime_zz,hprimewgll_zz,wxgll,wzgll, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right, &
- nspec_inner, ispec_inner_to_glob, .false. &
+ nspec_inner, ispec_inner_to_glob, .false. &
)
endif
@@ -2055,7 +2055,7 @@
endif
#endif
-! second call, computation on inner elements and update of
+! second call, computation on inner elements and update of
if(any_elastic) &
call compute_forces_elastic(npoin,nspec,nelemabs,numat,iglob_source, &
ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs,assign_external_model, &
@@ -2256,7 +2256,7 @@
endif
if(imagetype == 1) then
-
+
if ( myrank == 0 ) then
write(IOUT,*) 'drawing displacement vector as small arrows...'
endif
@@ -2362,7 +2362,7 @@
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
else if(imagetype == 3) then
-
+
if ( myrank == 0 ) then
write(IOUT,*) 'drawing image of vertical component of acceleration vector...'
endif
@@ -2386,7 +2386,7 @@
endif
image_color_data(:,:) = 0.d0
-
+
do k = 1, nb_pixel_loc
j = ceiling(real(num_pixel_loc(k)) / real(NX_IMAGE_color))
i = num_pixel_loc(k) - (j-1)*NX_IMAGE_color
Modified: seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2007-09-27 15:27:32 UTC (rev 8592)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2007-12-07 23:59:18 UTC (rev 8593)
@@ -4,10 +4,10 @@
! S P E C F E M 2 D Version 5.2
! ------------------------------
!
-! Dimitri Komatitsch
+! Main authors: Dimitri Komatitsch, Nicolas Le Goff and Roland Martin
! University of Pau, France
!
-! (c) April 2007
+! (c) November 2007
!
!========================================================================
@@ -26,10 +26,10 @@
#endif
integer :: nrec,NSTEP,it,seismotype
- integer :: NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current
+ integer :: NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current
double precision :: t0,deltat
-
+
integer, intent(in) :: nrecloc,myrank
integer, dimension(nrec),intent(in) :: which_proc_receiver
@@ -52,7 +52,7 @@
! scaling factor for Seismic Unix xsu dislay
double precision, parameter :: FACTORXSU = 1.d0
-
+
integer :: irecloc
#ifdef USE_MPI
@@ -87,28 +87,28 @@
allocate(buffer_binary(NTSTEP_BETWEEN_OUTPUT_SEISMO,number_of_components))
-
+
if ( myrank == 0 .and. seismo_offset == 0 ) then
-
+
! delete the old files
open(unit=11,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown')
close(11,status='delete')
-
+
open(unit=11,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown')
close(11,status='delete')
-
+
open(unit=11,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown')
close(11,status='delete')
-
+
open(unit=11,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown')
close(11,status='delete')
-
+
open(unit=11,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown')
close(11,status='delete')
open(unit=11,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown')
close(11,status='delete')
-
+
endif
if ( myrank == 0 ) then
@@ -130,25 +130,25 @@
if(seismotype /= 4) then
open(unit=14,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown',access='direct',recl=4)
open(unit=15,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown',access='direct',recl=8)
-
+
end if
-
+
end if
irecloc = 0
do irec = 1,nrec
-
+
if ( myrank == 0 ) then
-
+
if ( which_proc_receiver(irec) == myrank ) then
irecloc = irecloc + 1
buffer_binary(:,1) = sisux(:,irecloc)
if ( number_of_components == 2 ) then
buffer_binary(:,2) = sisuz(:,irecloc)
end if
-
-#ifdef USE_MPI
+
+#ifdef USE_MPI
else
call MPI_RECV(buffer_binary(1,1),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,&
which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
@@ -156,14 +156,14 @@
call MPI_RECV(buffer_binary(1,2),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,&
which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
end if
-
-
+
+
#endif
end if
-
+
! write trace
do iorientation = 1,number_of_components
-
+
if(iorientation == 1) then
chn = 'BHX'
else if(iorientation == 2) then
@@ -171,10 +171,10 @@
else
call exit_MPI('incorrect channel value')
endif
-
+
! in case of pressure, use different abbreviation
if(seismotype == 4) chn = 'PRE'
-
+
! create the name of the seismogram file for each slice
! file name includes the name of the station, the network and the component
length_station_name = len_trim(station_name(irec))
@@ -183,14 +183,14 @@
! check that length conforms to standard
if(length_station_name < 1 .or. length_station_name > MAX_LENGTH_STATION_NAME) then
call exit_MPI('wrong length of station name')
- end if
- if(length_network_name < 1 .or. length_network_name > MAX_LENGTH_NETWORK_NAME) then
+ end if
+ if(length_network_name < 1 .or. length_network_name > MAX_LENGTH_NETWORK_NAME) then
call exit_MPI('wrong length of network name')
end if
-
+
write(sisname,"('OUTPUT_FILES/',a,'.',a,'.',a3,'.sem',a1)") station_name(irec)(1:length_station_name),&
network_name(irec)(1:length_network_name),chn,component
-
+
! save seismograms in text format with no subsampling.
! Because we do not subsample the output, this can result in large files
! if the simulation uses many time steps. However, subsampling the output
@@ -201,7 +201,7 @@
close(11,status='delete')
endif
open(unit=11,file=sisname(1:len_trim(sisname)),status='unknown',position='append')
-
+
! make sure we never write more than the maximum number of time steps
! subtract offset of the source to make sure travel time is correct
do isample = 1,seismo_current
@@ -211,7 +211,7 @@
write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ',sngl(buffer_binary(isample,iorientation))
endif
enddo
-
+
close(11)
end do
@@ -232,11 +232,11 @@
call MPI_SEND(sisux(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
if ( number_of_components == 2 ) then
call MPI_SEND(sisuz(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
- end if
+ end if
end if
#endif
-
+
end if
enddo
@@ -254,7 +254,7 @@
!----
if ( myrank == 0 ) then
-
+
! ligne de recepteurs pour Xsu
open(unit=11,file='OUTPUT_FILES/receiver_line_Xsu_XWindow',status='unknown')
More information about the cig-commits
mailing list