[cig-commits] r18905 - in seismo/3D/SPECFEM3D/trunk/src: shared specfem3D
yangl at geodynamics.org
yangl at geodynamics.org
Wed Sep 14 13:28:51 PDT 2011
Author: yangl
Date: 2011-09-14 13:28:50 -0700 (Wed, 14 Sep 2011)
New Revision: 18905
Modified:
seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90
seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
seismo/3D/SPECFEM3D/trunk/src/shared/parallel.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90
Log:
(SPECFEM3D_SESAME) Seismic Unix format related, NOT benchmarked yet
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90 2011-09-14 18:52:37 UTC (rev 18904)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90 2011-09-14 20:28:50 UTC (rev 18905)
@@ -52,7 +52,7 @@
integer :: NSPEC_AB, NGLOB_AB
integer :: numpoin
integer :: i, ios, it, ier
- integer :: iproc, proc1, proc2, num_node, node_list(300)
+ integer :: iproc, proc1, proc2, num_node, node_list(2000)
integer :: np, ne, npp, nee, nelement, njunk
character(len=256) :: sline, arg(6), filename, indir, outdir
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2011-09-14 18:52:37 UTC (rev 18904)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/constants.h.in 2011-09-14 20:28:50 UTC (rev 18905)
@@ -75,6 +75,8 @@
! ouput format of seismograms, ASCII or binary
logical, parameter :: SEISMOGRAMS_BINARY = .false.
+! output format of seismograms, Seismic Unix (binary with 240-byte-headers)
+ logical, parameter :: SU_FORMAT=.false.
! input, output and main MPI I/O files
integer, parameter :: ISTANDARD_OUTPUT = 6
@@ -94,6 +96,10 @@
integer, parameter :: IIN_INTERFACES = 43
! I/O unit for noise simulations
integer, parameter :: IIN_NOISE = 44,IOUT_NOISE = 45
+! I/O unit for SU formatted seismograms
+ integer, parameter :: IOUT_SU = 46
+! I/O unit for SU formatted adjoint sources
+ integer, parameter :: IIN_SU1 = 47, IIN_SU2 = 48, IIN_SU3 = 49
! ignore variable name field (junk) at the beginning of each input line
logical, parameter :: IGNORE_JUNK = .true.,DONT_IGNORE_JUNK = .false.
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/parallel.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/parallel.f90 2011-09-14 18:52:37 UTC (rev 18904)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/parallel.f90 2011-09-14 20:28:50 UTC (rev 18905)
@@ -976,7 +976,47 @@
!----
!
+ subroutine send_dp(sendbuf, sendcount, dest, sendtag)
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ integer dest,sendtag
+ integer sendcount
+ double precision,dimension(sendcount):: sendbuf
+ integer ier
+
+ call MPI_SEND(sendbuf,sendcount,MPI_DOUBLE_PRECISION,dest,sendtag,MPI_COMM_WORLD,ier)
+
+ end subroutine send_dp
+
+!
+!----
+!
+
+ subroutine recv_dp(recvbuf, recvcount, dest, recvtag)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ integer dest,recvtag
+ integer recvcount
+ double precision,dimension(recvcount):: recvbuf
+ integer req(MPI_STATUS_SIZE)
+ integer ier
+
+ call MPI_RECV(recvbuf,recvcount,MPI_DOUBLE_PRECISION,dest,recvtag,MPI_COMM_WORLD,req,ier)
+
+ end subroutine recv_dp
+
+!
+!----
+!
+
subroutine sendv_cr(sendbuf, sendcount, dest, sendtag)
implicit none
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2011-09-14 18:52:37 UTC (rev 18904)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2011-09-14 20:28:50 UTC (rev 18905)
@@ -39,7 +39,7 @@
use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
xigll,yigll,zigll,xi_receiver,eta_receiver,gamma_receiver,&
- station_name,network_name,adj_source_file
+ station_name,network_name,adj_source_file,nrec_local,number_receiver_global
implicit none
include "constants.h"
@@ -89,6 +89,18 @@
integer :: isource,iglob,ispec,i,j,k,ier
integer :: irec_local,irec
+! adjoint sources in SU format
+ integer :: it_start,it_end
+ real(kind=CUSTOM_REAL) :: adj_temp(NSTEP)
+ real(kind=CUSTOM_REAL) :: adj_src(NTSTEP_BETWEEN_READ_ADJSRC,NDIM)
+ character(len=256) :: procname
+ integer,parameter :: nheader=240 ! 240 bytes
+ integer(kind=2) :: i2head(nheader/2) ! 2-byte-integer
+ integer(kind=4) :: i4head(nheader/4) ! 4-byte-integer
+ real(kind=4) :: r4head(nheader/4) ! 4-byte-real
+ equivalence (i2head,i4head,r4head) ! share the same 240-byte-memory
+ double precision :: hxir(NGLLX), hpxir(NGLLX), hetar(NGLLY), hpetar(NGLLY),hgammar(NGLLZ), hpgammar(NGLLZ)
+
! plotting source time function
if(PRINT_SOURCE_TIME_FUNCTION .and. .not. phase_is_inner ) then
! initializes total
@@ -235,24 +247,59 @@
allocate(adj_sourcearray(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
if( ier /= 0 ) stop 'error allocating array adj_sourcearray'
- irec_local = 0
- do irec = 1, nrec
- ! compute source arrays
- if (myrank == islice_selected_rec(irec)) then
- irec_local = irec_local + 1
+ if (.not. SU_FORMAT) then
+ !!! read ascii adjoint sources
+ irec_local = 0
+ do irec = 1, nrec
+ ! compute source arrays
+ if (myrank == islice_selected_rec(irec)) then
+ irec_local = irec_local + 1
- ! reads in **sta**.**net**.**LH**.adj files
- adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
- call compute_arrays_adjoint_source(myrank,adj_source_file, &
- xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
- adj_sourcearray, xigll,yigll,zigll, &
- it_sub_adj,NSTEP,NTSTEP_BETWEEN_READ_ADJSRC)
- do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
- adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
- enddo
+ ! reads in **sta**.**net**.**LH**.adj files
+ adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
+ call compute_arrays_adjoint_source(myrank,adj_source_file, &
+ xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
+ adj_sourcearray, xigll,yigll,zigll, &
+ it_sub_adj,NSTEP,NTSTEP_BETWEEN_READ_ADJSRC)
+ do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
+ adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
+ enddo
- endif
- enddo
+ endif
+ enddo
+ else
+ !!! read SU adjoint sources
+ ! range of the block we need to read
+ it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC + 1
+ it_end = it_start + NTSTEP_BETWEEN_READ_ADJSRC - 1
+ write(procname,"(i4)") myrank
+ open(unit=IIN_SU1, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj', &
+ access='direct',recl=240+4*(NSTEP))
+ do irec_local = 1,nrec_local
+ irec = number_receiver_global(irec_local)
+ read(IIN_SU1,rec=irec_local) r4head, adj_temp
+ adj_src(:,1)=adj_temp(it_start:it_end)
+ adj_src(:,2)=0.0 !TRIVIAL
+ adj_src(:,3)=0.0 !TRIVIAL
+ ! lagrange interpolators for receiver location
+ call lagrange_any(xi_receiver(irec),NGLLX,xigll,hxir,hpxir)
+ call lagrange_any(eta_receiver(irec),NGLLY,yigll,hetar,hpetar)
+ call lagrange_any(gamma_receiver(irec),NGLLZ,zigll,hgammar,hpgammar)
+ ! interpolates adjoint source onto GLL points within this element
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ adj_sourcearray(:,:,i,j,k) = hxir(i) * hetar(j) * hgammar(k) * adj_src(:,:)
+ enddo
+ enddo
+ enddo
+ do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
+ adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
+ enddo
+ enddo
+ close(IIN_SU1)
+ endif !if (.not. SU_FORMAT)
+
deallocate(adj_sourcearray)
endif ! if(ibool_read_adj_arrays)
@@ -276,13 +323,10 @@
do i = 1,NGLLX
iglob = ibool(i,j,k,ispec)
- ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
- ! note: it takes the first component of the adj_sourcearrays
- ! the idea is to have e.g. a pressure source, where all 3 components would be the same
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - adj_sourcearrays(irec_local, &
+ + adj_sourcearrays(irec_local, &
NTSTEP_BETWEEN_READ_ADJSRC - mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC), &
- 1,i,j,k) / kappastore(i,j,k,ispec)
+ 1,i,j,k)
enddo
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2011-09-14 18:52:37 UTC (rev 18904)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2011-09-14 20:28:50 UTC (rev 18905)
@@ -42,7 +42,8 @@
LOCAL_PATH,wgllwgll_xy, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk,free_surface_jacobian2Dw, &
noise_sourcearray,irec_master_noise, &
- normal_x_noise,normal_y_noise,normal_z_noise,mask_noise,noise_surface_movie
+ normal_x_noise,normal_y_noise,normal_z_noise,mask_noise,noise_surface_movie, &
+ nrec_local,number_receiver_global
use specfem_par_movie,only: &
store_val_ux_external_mesh,store_val_uy_external_mesh,store_val_uz_external_mesh
@@ -94,6 +95,18 @@
integer :: isource,iglob,i,j,k,ispec,ier
integer :: irec_local,irec
+! adjoint sources in SU format
+ integer :: it_start,it_end
+ real(kind=CUSTOM_REAL) :: adj_temp(NSTEP)
+ real(kind=CUSTOM_REAL) :: adj_src(NTSTEP_BETWEEN_READ_ADJSRC,NDIM)
+ character(len=256) :: procname
+ integer,parameter :: nheader=240 ! 240 bytes
+ integer(kind=2) :: i2head(nheader/2) ! 2-byte-integer
+ integer(kind=4) :: i4head(nheader/4) ! 4-byte-integer
+ real(kind=4) :: r4head(nheader/4) ! 4-byte-real
+ equivalence (i2head,i4head,r4head) ! share the same 240-byte-memory
+ double precision :: hxir(NGLLX), hpxir(NGLLX), hetar(NGLLY), hpetar(NGLLY),hgammar(NGLLZ), hpgammar(NGLLZ)
+
! plotting source time function
if(PRINT_SOURCE_TIME_FUNCTION .and. .not. phase_is_inner ) then
! initializes total
@@ -218,25 +231,68 @@
! allocates temporary source array
allocate(adj_sourcearray(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
if( ier /= 0 ) stop 'error allocating array adj_sourcearray'
-
- irec_local = 0
- do irec = 1, nrec
- ! compute source arrays
- if (myrank == islice_selected_rec(irec)) then
- irec_local = irec_local + 1
- ! reads in **sta**.**net**.**LH**.adj files
- adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
- call compute_arrays_adjoint_source(myrank,adj_source_file, &
- xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
- adj_sourcearray, xigll,yigll,zigll, &
- it_sub_adj,NSTEP,NTSTEP_BETWEEN_READ_ADJSRC)
- do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
- adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
- enddo
+ if (.not. SU_FORMAT) then
+ !!! read ascii adjoint sources
+ irec_local = 0
+ do irec = 1, nrec
+ ! compute source arrays
+ if (myrank == islice_selected_rec(irec)) then
+ irec_local = irec_local + 1
+ ! reads in **sta**.**net**.**LH**.adj files
+ adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
+ call compute_arrays_adjoint_source(myrank,adj_source_file, &
+ xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
+ adj_sourcearray, xigll,yigll,zigll, &
+ it_sub_adj,NSTEP,NTSTEP_BETWEEN_READ_ADJSRC)
+ do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
+ adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
+ enddo
- endif
- enddo
+ endif
+ enddo
+ else
+ !!! read SU adjoint sources
+ ! range of the block we need to read
+ it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC + 1
+ it_end = it_start + NTSTEP_BETWEEN_READ_ADJSRC - 1
+ write(procname,"(i4)") myrank
+ ! read adjoint sources
+ open(unit=IIN_SU1, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dx_SU.adj', &
+ access='direct',recl=240+4*NSTEP)
+ open(unit=IIN_SU2, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dy_SU.adj', &
+ access='direct',recl=240+4*NSTEP)
+ open(unit=IIN_SU3, file=trim(adjustl(OUTPUT_FILES_PATH))//'../SEM/'//trim(adjustl(procname))//'_dz_SU.adj', &
+ access='direct',recl=240+4*NSTEP)
+ do irec_local = 1,nrec_local
+ irec = number_receiver_global(irec_local)
+ read(IIN_SU1,rec=irec_local) r4head, adj_temp
+ adj_src(:,1)=adj_temp(it_start:it_end)
+ read(IIN_SU2,rec=irec_local) r4head, adj_temp
+ adj_src(:,2)=adj_temp(it_start:it_end)
+ read(IIN_SU3,rec=irec_local) r4head, adj_temp
+ adj_src(:,3)=adj_temp(it_start:it_end)
+ ! lagrange interpolators for receiver location
+ call lagrange_any(xi_receiver(irec),NGLLX,xigll,hxir,hpxir)
+ call lagrange_any(eta_receiver(irec),NGLLY,yigll,hetar,hpetar)
+ call lagrange_any(gamma_receiver(irec),NGLLZ,zigll,hgammar,hpgammar)
+ ! interpolates adjoint source onto GLL points within this element
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ adj_sourcearray(:,:,i,j,k) = hxir(i) * hetar(j) * hgammar(k) * adj_src(:,:)
+ enddo
+ enddo
+ enddo
+ do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
+ adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
+ enddo
+ enddo
+ close(IIN_SU1)
+ close(IIN_SU2)
+ close(IIN_SU3)
+ endif !if(.not. SU_FORMAT)
+
deallocate(adj_sourcearray)
endif ! if(ibool_read_adj_arrays)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2011-09-14 18:52:37 UTC (rev 18904)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2011-09-14 20:28:50 UTC (rev 18905)
@@ -73,7 +73,6 @@
double precision, allocatable, dimension(:) :: x_target,y_target,z_target
double precision, allocatable, dimension(:) :: horiz_dist
double precision, allocatable, dimension(:) :: x_found,y_found,z_found
- double precision, allocatable, dimension(:,:) :: x_found_all,y_found_all,z_found_all
integer irec
integer i,j,k,ispec,iglob,iface,inode,imin,imax,jmin,jmax,kmin,kmax,igll,jgll,kgll
@@ -113,7 +112,6 @@
! use dynamic allocation
double precision, dimension(:), allocatable :: final_distance
- double precision, dimension(:,:), allocatable :: final_distance_all
double precision distmin,final_distance_max
! receiver information
@@ -126,17 +124,37 @@
double precision, dimension(3,3,nrec) :: nu
character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name
character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name
-
- integer, allocatable, dimension(:,:) :: ispec_selected_rec_all
double precision, allocatable, dimension(:) :: stlat,stlon,stele,stbur,stutm_x,stutm_y,elevation
- double precision, allocatable, dimension(:,:) :: xi_receiver_all,eta_receiver_all,gamma_receiver_all
- double precision, allocatable, dimension(:,:,:,:) :: nu_all
+ double precision, allocatable, dimension(:) :: x_found_all,y_found_all,z_found_all
+ double precision, dimension(:), allocatable :: final_distance_all
+ integer, allocatable, dimension(:) :: ispec_selected_rec_all
+ double precision, allocatable, dimension(:) :: xi_receiver_all,eta_receiver_all,gamma_receiver_all
+ double precision, allocatable, dimension(:,:,:) :: nu_all
+
integer :: ier
character(len=256) OUTPUT_FILES
+ real(kind=CUSTOM_REAL) :: xmin,xmax,ymin,ymax,zmin,zmax
+ real(kind=CUSTOM_REAL) :: xmin_ELE,xmax_ELE,ymin_ELE,ymax_ELE,zmin_ELE,zmax_ELE
+ integer :: imin_temp,imax_temp,jmin_temp,jmax_temp,kmin_temp,kmax_temp
+ integer,dimension(NGLLX*NGLLY*NGLLZ) :: iglob_temp
+
! **************
+ ! dimension of model in current proc
+ xmin=minval(xstore(:)); xmax=maxval(xstore(:))
+ ymin=minval(ystore(:)); ymax=maxval(ystore(:))
+ zmin=minval(zstore(:)); zmax=maxval(zstore(:))
+ if(FASTER_RECEIVERS_POINTS_ONLY) then
+ imin_temp = 1; imax_temp = NGLLX
+ jmin_temp = 1; jmax_temp = NGLLY
+ kmin_temp = 1; kmax_temp = NGLLZ
+ else
+ imin_temp = 2; imax_temp = NGLLX - 1
+ jmin_temp = 2; jmax_temp = NGLLY - 1
+ kmin_temp = 2; kmax_temp = NGLLZ - 1
+ endif
! get MPI starting time
time_start = wtime()
@@ -182,15 +200,15 @@
y_found(nrec), &
z_found(nrec), &
final_distance(nrec), &
- ispec_selected_rec_all(nrec,0:NPROC-1), &
- xi_receiver_all(nrec,0:NPROC-1), &
- eta_receiver_all(nrec,0:NPROC-1), &
- gamma_receiver_all(nrec,0:NPROC-1), &
- x_found_all(nrec,0:NPROC-1), &
- y_found_all(nrec,0:NPROC-1), &
- z_found_all(nrec,0:NPROC-1), &
- final_distance_all(nrec,0:NPROC-1), &
- nu_all(3,3,nrec,0:NPROC-1),stat=ier)
+ ispec_selected_rec_all(nrec), &
+ xi_receiver_all(nrec), &
+ eta_receiver_all(nrec), &
+ gamma_receiver_all(nrec), &
+ x_found_all(nrec), &
+ y_found_all(nrec), &
+ z_found_all(nrec), &
+ final_distance_all(nrec), &
+ nu_all(3,3,nrec),stat=ier)
if( ier /= 0 ) stop 'error allocating arrays for locating receivers'
! loop on all the stations
@@ -325,78 +343,136 @@
endif
!if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
- ! determines closest GLL point
- ispec_selected_rec(irec) = 0
- do ispec=1,NSPEC_AB
+ if (.not. SU_FORMAT) then
+ ! determines closest GLL point
+ ispec_selected_rec(irec) = 0
+ do ispec=1,NSPEC_AB
- ! define the interval in which we look for points
- if(FASTER_RECEIVERS_POINTS_ONLY) then
- imin = 1
- imax = NGLLX
+ ! define the interval in which we look for points
+ if(FASTER_RECEIVERS_POINTS_ONLY) then
+ imin = 1
+ imax = NGLLX
- jmin = 1
- jmax = NGLLY
+ jmin = 1
+ jmax = NGLLY
- kmin = 1
- kmax = NGLLZ
+ kmin = 1
+ kmax = NGLLZ
- else
- ! loop only on points inside the element
- ! exclude edges to ensure this point is not shared with other elements
- imin = 2
- imax = NGLLX - 1
+ else
+ ! loop only on points inside the element
+ ! exclude edges to ensure this point is not shared with other elements
+ imin = 2
+ imax = NGLLX - 1
- jmin = 2
- jmax = NGLLY - 1
+ jmin = 2
+ jmax = NGLLY - 1
- kmin = 2
- kmax = NGLLZ - 1
- endif
+ kmin = 2
+ kmax = NGLLZ - 1
+ endif
- do k = kmin,kmax
- do j = jmin,jmax
- do i = imin,imax
+ do k = kmin,kmax
+ do j = jmin,jmax
+ do i = imin,imax
- iglob = ibool(i,j,k,ispec)
+ iglob = ibool(i,j,k,ispec)
- if (.not. RECVS_CAN_BE_BURIED_EXT_MESH) then
- if ((.not. iglob_is_surface_external_mesh(iglob)) .or. (.not. ispec_is_surface_external_mesh(ispec))) then
- cycle
- endif
- endif
+ if (.not. RECVS_CAN_BE_BURIED_EXT_MESH) then
+ if ((.not. iglob_is_surface_external_mesh(iglob)) .or. (.not. ispec_is_surface_external_mesh(ispec))) then
+ cycle
+ endif
+ endif
- dist = dsqrt((x_target(irec)-dble(xstore(iglob)))**2 &
- +(y_target(irec)-dble(ystore(iglob)))**2 &
- +(z_target(irec)-dble(zstore(iglob)))**2)
+ dist = dsqrt((x_target(irec)-dble(xstore(iglob)))**2 &
+ +(y_target(irec)-dble(ystore(iglob)))**2 &
+ +(z_target(irec)-dble(zstore(iglob)))**2)
- ! keep this point if it is closer to the receiver
- if(dist < distmin) then
- distmin = dist
- ispec_selected_rec(irec) = ispec
- ix_initial_guess(irec) = i
- iy_initial_guess(irec) = j
- iz_initial_guess(irec) = k
+ ! keep this point if it is closer to the receiver
+ if(dist < distmin) then
+ distmin = dist
+ ispec_selected_rec(irec) = ispec
+ ix_initial_guess(irec) = i
+ iy_initial_guess(irec) = j
+ iz_initial_guess(irec) = k
- xi_receiver(irec) = dble(ix_initial_guess(irec))
- eta_receiver(irec) = dble(iy_initial_guess(irec))
- gamma_receiver(irec) = dble(iz_initial_guess(irec))
- x_found(irec) = xstore(iglob)
- y_found(irec) = ystore(iglob)
- z_found(irec) = zstore(iglob)
- endif
+ xi_receiver(irec) = dble(ix_initial_guess(irec))
+ eta_receiver(irec) = dble(iy_initial_guess(irec))
+ gamma_receiver(irec) = dble(iz_initial_guess(irec))
+ x_found(irec) = xstore(iglob)
+ y_found(irec) = ystore(iglob)
+ z_found(irec) = zstore(iglob)
+ endif
- enddo
- enddo
- enddo
+ enddo
+ enddo
+ enddo
- ! compute final distance between asked and found (converted to km)
- final_distance(irec) = dsqrt((x_target(irec)-x_found(irec))**2 + &
- (y_target(irec)-y_found(irec))**2 + (z_target(irec)-z_found(irec))**2)
+ ! compute final distance between asked and found (converted to km)
+ final_distance(irec) = dsqrt((x_target(irec)-x_found(irec))**2 + &
+ (y_target(irec)-y_found(irec))**2 + (z_target(irec)-z_found(irec))**2)
- ! end of loop on all the spectral elements in current slice
- enddo
+ ! end of loop on all the spectral elements in current slice
+ enddo
+ else
+ ispec_selected_rec(irec) = 0
+ ix_initial_guess(irec) = 0
+ iy_initial_guess(irec) = 0
+ iz_initial_guess(irec) = 0
+ final_distance(irec) = HUGEVAL
+ if ( (x_target(irec)>=xmin .and. x_target(irec)<=xmax) .and. &
+ (y_target(irec)>=ymin .and. y_target(irec)<=ymax) .and. &
+ (z_target(irec)>=zmin .and. z_target(irec)<=zmax) ) then
+ do ispec=1,NSPEC_AB
+ iglob_temp=reshape(ibool(:,:,:,ispec),(/NGLLX*NGLLY*NGLLZ/))
+ xmin_ELE=minval(xstore(iglob_temp))
+ xmax_ELE=maxval(xstore(iglob_temp))
+ ymin_ELE=minval(ystore(iglob_temp))
+ ymax_ELE=maxval(ystore(iglob_temp))
+ zmin_ELE=minval(zstore(iglob_temp))
+ zmax_ELE=maxval(zstore(iglob_temp))
+ if ( (x_target(irec)>=xmin_ELE .and. x_target(irec)<=xmax_ELE) .and. &
+ (y_target(irec)>=ymin_ELE .and. y_target(irec)<=ymax_ELE) .and. &
+ (z_target(irec)>=zmin_ELE .and. z_target(irec)<=zmax_ELE) ) then
+ ! we find the element (ispec) which "may" contain the receiver (irec)
+ ! so we only need to compute distances (which is expensive because of "dsqrt") within those elements
+ ispec_selected_rec(irec) = ispec
+ do k = kmin_temp,kmax_temp
+ do j = jmin_temp,jmax_temp
+ do i = imin_temp,imax_temp
+ iglob = ibool(i,j,k,ispec)
+ ! for comparison purpose, we don't have to do "dsqrt", which is expensive
+ dist = ((x_target(irec)-dble(xstore(iglob)))**2 &
+ +(y_target(irec)-dble(ystore(iglob)))**2 &
+ +(z_target(irec)-dble(zstore(iglob)))**2)
+ if(dist < distmin) then
+ distmin = dist
+ ix_initial_guess(irec) = i
+ iy_initial_guess(irec) = j
+ iz_initial_guess(irec) = k
+ xi_receiver(irec) = dble(ix_initial_guess(irec))
+ eta_receiver(irec) = dble(iy_initial_guess(irec))
+ gamma_receiver(irec) = dble(iz_initial_guess(irec))
+ x_found(irec) = xstore(iglob)
+ y_found(irec) = ystore(iglob)
+ z_found(irec) = zstore(iglob)
+ endif
+ enddo
+ enddo
+ enddo
+ final_distance(irec) = dsqrt((x_target(irec)-x_found(irec))**2 + &
+ (y_target(irec)-y_found(irec))**2 + (z_target(irec)-z_found(irec))**2)
+ endif ! if receiver "may" be within this element
+ enddo ! do ispec=1,NSPEC_AB
+ endif ! if receiver "may" be within this proc
+ endif !if (.not. SU_FORMAT)
if (ispec_selected_rec(irec) == 0) then
+ ! receiver is NOT within this proc, assign trivial values
+ ispec_selected_rec(irec) = 1
+ ix_initial_guess(irec) = 1
+ iy_initial_guess(irec) = 1
+ iz_initial_guess(irec) = 1
final_distance(irec) = HUGEVAL
endif
@@ -684,46 +760,50 @@
! synchronize all the processes to make sure all the estimates are available
call sync_all()
-
! for MPI version, gather information from all the nodes
- ispec_selected_rec_all(:,:) = -1
- call gather_all_i(ispec_selected_rec,nrec,ispec_selected_rec_all,nrec,NPROC)
- call gather_all_dp(xi_receiver,nrec,xi_receiver_all,nrec,NPROC)
- call gather_all_dp(eta_receiver,nrec,eta_receiver_all,nrec,NPROC)
- call gather_all_dp(gamma_receiver,nrec,gamma_receiver_all,nrec,NPROC)
- call gather_all_dp(final_distance,nrec,final_distance_all,nrec,NPROC)
- call gather_all_dp(x_found,nrec,x_found_all,nrec,NPROC)
- call gather_all_dp(y_found,nrec,y_found_all,nrec,NPROC)
- call gather_all_dp(z_found,nrec,z_found_all,nrec,NPROC)
- call gather_all_dp(nu,3*3*nrec,nu_all,3*3*nrec,NPROC)
+ if (myrank/=0) then ! gather information from other processors (one at a time)
+ call send_i(ispec_selected_rec, nrec,0,0)
+ call send_dp(xi_receiver, nrec,0,1)
+ call send_dp(eta_receiver, nrec,0,2)
+ call send_dp(gamma_receiver, nrec,0,3)
+ call send_dp(final_distance, nrec,0,4)
+ call send_dp(x_found, nrec,0,5)
+ call send_dp(y_found, nrec,0,6)
+ call send_dp(z_found, nrec,0,7)
+ call send_dp(nu, 3*3*nrec,0,8)
+ else
+ islice_selected_rec(:) = 0
+ do iprocloop=1,NPROC-1
+ call recv_i(ispec_selected_rec_all, nrec,iprocloop,0)
+ call recv_dp(xi_receiver_all, nrec,iprocloop,1)
+ call recv_dp(eta_receiver_all, nrec,iprocloop,2)
+ call recv_dp(gamma_receiver_all, nrec,iprocloop,3)
+ call recv_dp(final_distance_all, nrec,iprocloop,4)
+ call recv_dp(x_found_all, nrec,iprocloop,5)
+ call recv_dp(y_found_all, nrec,iprocloop,6)
+ call recv_dp(z_found_all, nrec,iprocloop,7)
+ call recv_dp(nu_all, 3*3*nrec,iprocloop,8)
+ do irec=1,nrec
+ if (final_distance_all(irec) < final_distance(irec)) then
+ final_distance(irec) = final_distance_all(irec)
+ islice_selected_rec(irec) = iprocloop
+ ispec_selected_rec(irec) = ispec_selected_rec_all(irec)
+ xi_receiver(irec) = xi_receiver_all(irec)
+ eta_receiver(irec) = eta_receiver_all(irec)
+ gamma_receiver(irec) = gamma_receiver_all(irec)
+ x_found(irec) = x_found_all(irec)
+ y_found(irec) = y_found_all(irec)
+ z_found(irec) = z_found_all(irec)
+ nu(:,:,irec) = nu_all(:,:,irec)
+ endif
+ enddo
+ enddo
+ endif
+ call sync_all()
! this is executed by main process only
if(myrank == 0) then
- ! check that the gather operation went well
- if(any(ispec_selected_rec_all(:,:) == -1)) call exit_MPI(myrank,'gather operation failed for receivers')
-
- ! MPI loop on all the results to determine the best slice
- islice_selected_rec(:) = -1
- do irec = 1,nrec
- distmin = HUGEVAL
- do iprocloop = 0,NPROC-1
- if(final_distance_all(irec,iprocloop) < distmin) then
- distmin = final_distance_all(irec,iprocloop)
- islice_selected_rec(irec) = iprocloop
- ispec_selected_rec(irec) = ispec_selected_rec_all(irec,iprocloop)
- xi_receiver(irec) = xi_receiver_all(irec,iprocloop)
- eta_receiver(irec) = eta_receiver_all(irec,iprocloop)
- gamma_receiver(irec) = gamma_receiver_all(irec,iprocloop)
- x_found(irec) = x_found_all(irec,iprocloop)
- y_found(irec) = y_found_all(irec,iprocloop)
- z_found(irec) = z_found_all(irec,iprocloop)
- nu(:,:,irec) = nu_all(:,:,irec,iprocloop)
- endif
- enddo
- final_distance(irec) = distmin
- enddo
-
do irec=1,nrec
write(IMAIN,*)
@@ -809,12 +889,19 @@
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', OUTPUT_FILES_PATH(1:len_trim(OUTPUT_FILES_PATH)))
- ! write the list of stations and associated epicentral distance
- open(unit=27,file=trim(OUTPUT_FILES)//'/output_list_stations.txt',status='unknown')
+ !! write the list of stations and associated epicentral distance
+ !open(unit=27,file=trim(OUTPUT_FILES)//'/output_list_stations.txt',status='unknown')
+ !do irec=1,nrec
+ ! write(27,*) station_name(irec),'.',network_name(irec),' : ',horiz_dist(irec),' km horizontal distance'
+ !enddo
+ !close(27)
+
+ ! write the locations of stations, so that we can load them and write them to SU headers later
+ open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/output_list_stations.txt',status='unknown')
do irec=1,nrec
- write(27,*) station_name(irec),'.',network_name(irec),' : ',horiz_dist(irec),' km horizontal distance'
+ write(IOUT_SU,*) x_found(irec),y_found(irec),z_found(irec)
enddo
- close(27)
+ close(IOUT_SU)
! elapsed time since beginning of mesh generation
tCPU = wtime() - time_start
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2011-09-14 18:52:37 UTC (rev 18904)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2011-09-14 20:28:50 UTC (rev 18905)
@@ -936,6 +936,12 @@
write(IMAIN,*)
write(IMAIN,*) 'End of source detection - done'
write(IMAIN,*)
+ ! output source information to a file so that we can load it and write to SU headers later
+ open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/output_list_sources.txt',status='unknown')
+ do isource=1,NSOURCES
+ write(IOUT_SU,*) x_found_source(isource),y_found_source(isource),z_found_source(isource)
+ enddo
+ close(IOUT_SU)
endif
end subroutine locate_source
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2011-09-14 18:52:37 UTC (rev 18904)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2011-09-14 20:28:50 UTC (rev 18905)
@@ -571,7 +571,7 @@
! ADJOINT simulations
if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
-
+ if (.not.SU_FORMAT) then
! gets channel names
do icomp=1,NDIM
call write_channel_name(icomp,comp(icomp))
@@ -631,6 +631,14 @@
allocate(adj_sourcearrays(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
if( ier /= 0 ) stop 'error allocating array adj_sourcearrays'
adj_sourcearrays = 0._CUSTOM_REAL
+ else
+ ! skip counting, because only one file per component per proc in SU_FORMAT
+ nadj_rec_local=nrec_local
+ nadj_files_found=nrec_local
+ allocate(adj_sourcearrays(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array adj_sourcearrays'
+ adj_sourcearrays = 0._CUSTOM_REAL
+ endif !if (.not. SU_FORMAT)
else
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90 2011-09-14 18:52:37 UTC (rev 18904)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90 2011-09-14 20:28:50 UTC (rev 18905)
@@ -233,7 +233,7 @@
enddo ! nrec_local
! write the current or final seismograms
- if(mod(it,NTSTEP_BETWEEN_OUTPUT_SEISMOS) == 0 .or. it == NSTEP) then
+ if((mod(it,NTSTEP_BETWEEN_OUTPUT_SEISMOS) == 0 .or. it == NSTEP) .and. (.not.SU_FORMAT)) then
if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then
call write_seismograms_to_file(myrank,seismograms_d,number_receiver_global,station_name, &
network_name,nrec,nrec_local,it,DT,NSTEP,t0,LOCAL_PATH,1,SIMULATION_TYPE)
@@ -247,6 +247,11 @@
endif
endif
+! write ONE binary file for all receivers (nrec_local) within one proc
+! SU format, with 240-byte-header for each trace
+ if ((mod(it,NTSTEP_BETWEEN_OUTPUT_SEISMOS) == 0 .or. it==NSTEP) .and. SU_FORMAT) &
+ call write_seismograms_su()
+
end subroutine write_seismograms
@@ -747,3 +752,105 @@
endif
end subroutine band_instrument_code
+
+!=====================================================================
+
+ subroutine write_seismograms_su()
+
+ use specfem_par
+ use specfem_par_acoustic
+ use specfem_par_elastic
+ use specfem_par_poroelastic
+
+ implicit none
+
+ character(len=256) procname,final_LOCAL_PATH
+ integer :: irec_local,irec,ios
+ real :: x_station,y_station
+
+ ! headers
+ integer,parameter :: nheader=240 ! 240 bytes
+ integer(kind=2) :: i2head(nheader/2) ! 2-byte-integer
+ integer(kind=4) :: i4head(nheader/4) ! 4-byte-integer
+ real(kind=4) :: r4head(nheader/4) ! 4-byte-real
+ equivalence (i2head,i4head,r4head) ! share the same 240-byte-memory
+
+ double precision, allocatable, dimension(:) :: stlat,stlon,stele,stbur,stutm_x,stutm_y,elevation
+ double precision, allocatable, dimension(:) :: x_found,y_found,z_found
+ double precision :: x_found_source,y_found_source,z_found_source
+
+ allocate(x_found(nrec))
+ allocate(y_found(nrec))
+ allocate(z_found(nrec))
+ open(unit=IIN_SU1,file=trim(OUTPUT_FILES)//'/output_list_stations.txt',status='unknown')
+ do irec=1,nrec
+ read(IIN_SU1,*) x_found(irec),y_found(irec),z_found(irec)
+ enddo
+ close(IIN_SU1)
+ open(unit=IIN_SU1,file=trim(OUTPUT_FILES)//'/output_list_sources.txt',status='unknown')
+ read(IIN_SU1,*) x_found_source,y_found_source,z_found_source
+ close(IIN_SU1)
+ ! directory to store seismograms
+ if( USE_OUTPUT_FILES_PATH ) then
+ final_LOCAL_PATH = OUTPUT_FILES_PATH(1:len_trim(OUTPUT_FILES_PATH)) // '/'
+ else
+ ! create full final local path
+ final_LOCAL_PATH = trim(adjustl(LOCAL_PATH)) // '/'
+ endif
+ write(procname,"(i4)") myrank
+
+ ! write seismograms (dx)
+ open(unit=IOUT_SU, file=trim(adjustl(final_LOCAL_PATH))//trim(adjustl(procname))//'_dx_SU' ,&
+ form='unformatted', access='direct', recl=240+4*(NSTEP))
+ do irec_local = 1,nrec_local
+ irec = number_receiver_global(irec_local)
+ i4head(1) =irec
+ i4head(11) =z_found(irec)
+ i4head(13) =z_found_source
+ i4head(19) =x_found_source !utm_x_source(1)
+ i4head(20) =y_found_source !utm_y_source(1)
+ i4head(21) =x_found(irec) !stutm_x(irec)
+ i4head(22) =y_found(irec) !stutm_y(irec)
+ i2head(58) =NSTEP
+ i2head(59) =DT*1.0d6
+ write(IOUT_SU,rec=irec_local) r4head, seismograms_d(1,irec_local,:)
+ enddo
+ close(IOUT_SU)
+ ! write seismograms (dy)
+ open(unit=IOUT_SU, file=trim(adjustl(final_LOCAL_PATH))//trim(adjustl(procname))//'_dy_SU' ,&
+ form='unformatted', access='direct', recl=240+4*(NSTEP))
+ do irec_local = 1,nrec_local
+ irec = number_receiver_global(irec_local)
+ i4head(1) =irec
+ i4head(11) =z_found(irec)
+ i4head(13) =z_found_source
+ i4head(19) =x_found_source !utm_x_source(1)
+ i4head(20) =y_found_source !utm_y_source(1)
+ i4head(21) =x_found(irec) !stutm_x(irec)
+ i4head(22) =y_found(irec) !stutm_y(irec)
+ i2head(58) =NSTEP
+ i2head(59) =DT*1.0d6
+ write(IOUT_SU,rec=irec_local) r4head, seismograms_d(2,irec_local,:)
+ enddo
+ close(IOUT_SU)
+
+ ! write seismograms (dz)
+ open(unit=IOUT_SU, file=trim(adjustl(final_LOCAL_PATH))//trim(adjustl(procname))//'_dz_SU' ,&
+ form='unformatted', access='direct', recl=240+4*(NSTEP))
+ do irec_local = 1,nrec_local
+ irec = number_receiver_global(irec_local)
+ i4head(1) =irec
+ i4head(11) =z_found(irec)
+ i4head(13) =z_found_source
+ i4head(19) =x_found_source !utm_x_source(1)
+ i4head(20) =y_found_source !utm_y_source(1)
+ i4head(21) =x_found(irec) !stutm_x(irec)
+ i4head(22) =y_found(irec) !stutm_y(irec)
+ i2head(58) =NSTEP
+ i2head(59) =DT*1.0d6
+ write(IOUT_SU,rec=irec_local) r4head, seismograms_d(3,irec_local,:)
+ enddo
+ close(IOUT_SU)
+
+ end subroutine write_seismograms_su
+
More information about the CIG-COMMITS
mailing list