[cig-commits] r18878 - in seismo/2D/SPECFEM2D/trunk: DATA setup src/specfem2D
yangl at geodynamics.org
yangl at geodynamics.org
Wed Sep 7 10:23:21 PDT 2011
Author: yangl
Date: 2011-09-07 10:23:20 -0700 (Wed, 07 Sep 2011)
New Revision: 18878
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file.in
seismo/2D/SPECFEM2D/trunk/setup/constants.h.in
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_energy.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/paco_beyond_critical.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90
Log:
(NOT benchmarked yet) SPECFEM2D: seismograms in SeismicUnix format, acoustic kernels, inconsistent arguments fixes, minor bugs, etc.
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file.in 2011-09-07 17:01:35 UTC (rev 18877)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file.in 2011-09-07 17:23:20 UTC (rev 18878)
@@ -14,7 +14,7 @@
initialfield = .false. # use a plane wave as source or not
add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
assign_external_model = .false. # define external earth model or not
-READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+READ_EXTERNAL_SEP_FILE = .false. # Read external model from DATA/model_velocity.dat_input, or use routine
TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
TURN_VISCATTENUATION_ON = .false. # turn viscous attenuation on or off
Q0 = 1 # quality factor for viscous attenuation
@@ -34,7 +34,7 @@
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
# receiver line parameters for seismograms
-seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
+seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure 6=potential
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
Modified: seismo/2D/SPECFEM2D/trunk/setup/constants.h.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/setup/constants.h.in 2011-09-07 17:01:35 UTC (rev 18877)
+++ seismo/2D/SPECFEM2D/trunk/setup/constants.h.in 2011-09-07 17:23:20 UTC (rev 18878)
@@ -34,6 +34,9 @@
! this flag is ignored in the case of a serial simulation
logical, parameter :: FURTHER_REDUCE_CACHE_MISSES = .true.
+! output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+ logical,parameter :: SU_FORMAT=.false.
+
! for inverse Cuthill-McKee (1969) permutation
logical, parameter :: PERFORM_CUTHILL_MCKEE = .true.
logical, parameter :: INVERSE = .true.
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_energy.f90 2011-09-07 17:01:35 UTC (rev 18877)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_energy.f90 2011-09-07 17:23:20 UTC (rev 18878)
@@ -65,7 +65,7 @@
integer :: nspec,numat
! vector field in an element
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLX) :: vector_field_element
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLX) :: vector_field_element
! pressure in an element
integer :: N_SLS
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2011-09-07 17:01:35 UTC (rev 18877)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2011-09-07 17:23:20 UTC (rev 18878)
@@ -89,10 +89,10 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,rhoext
- double precision, dimension(NGLLZ,nspec_left,NSTEP) :: b_absorb_acoustic_left
- double precision, dimension(NGLLZ,nspec_right,NSTEP) :: b_absorb_acoustic_right
- double precision, dimension(NGLLX,nspec_top,NSTEP) :: b_absorb_acoustic_top
- double precision, dimension(NGLLX,nspec_bottom,NSTEP) :: b_absorb_acoustic_bottom
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,nspec_left,NSTEP) :: b_absorb_acoustic_left
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,nspec_right,NSTEP) :: b_absorb_acoustic_right
+ real(kind=CUSTOM_REAL), dimension(NGLLX,nspec_top,NSTEP) :: b_absorb_acoustic_top
+ real(kind=CUSTOM_REAL), dimension(NGLLX,nspec_bottom,NSTEP) :: b_absorb_acoustic_bottom
! derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
@@ -444,7 +444,7 @@
SIMULATION_TYPE,SAVE_FORWARD,nspec_left,nspec_right,&
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
b_absorb_acoustic_left,b_absorb_acoustic_right, &
- b_absorb_acoustic_bottom,b_absorb_acoustic_top)
+ b_absorb_acoustic_bottom,b_absorb_acoustic_top,IS_BACKWARD_FIELD)
! compute forces for the acoustic elements
@@ -471,7 +471,7 @@
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,rhoext
logical :: anyabs,assign_external_model
- logical :: SAVE_FORWARD
+ logical :: SAVE_FORWARD,IS_BACKWARD_FIELD
! derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
@@ -487,10 +487,10 @@
integer, dimension(nelemabs) :: ib_bottom
integer, dimension(nelemabs) :: ib_top
- double precision, dimension(NGLLZ,nspec_left,NSTEP) :: b_absorb_acoustic_left
- double precision, dimension(NGLLZ,nspec_right,NSTEP) :: b_absorb_acoustic_right
- double precision, dimension(NGLLX,nspec_top,NSTEP) :: b_absorb_acoustic_top
- double precision, dimension(NGLLX,nspec_bottom,NSTEP) :: b_absorb_acoustic_bottom
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,nspec_left,NSTEP) :: b_absorb_acoustic_left
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,nspec_right,NSTEP) :: b_absorb_acoustic_right
+ real(kind=CUSTOM_REAL), dimension(NGLLX,nspec_top,NSTEP) :: b_absorb_acoustic_top
+ real(kind=CUSTOM_REAL), dimension(NGLLX,nspec_bottom,NSTEP) :: b_absorb_acoustic_bottom
!---
!--- local variables
@@ -624,10 +624,17 @@
potential_dot_dot_acoustic(iglob) &
- potential_dot_acoustic(iglob)*weight/cpl/rhol
elseif(SIMULATION_TYPE == 2) then
- ! adds (previously) stored contribution
- potential_dot_dot_acoustic(iglob) = &
- potential_dot_dot_acoustic(iglob) &
- - b_absorb_acoustic_left(j,ib_left(ispecabs),NSTEP-it+1)
+ if(IS_BACKWARD_FIELD) then
+ ! adds (previously) stored contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - b_absorb_acoustic_left(j,ib_left(ispecabs),NSTEP-it+1)
+ else
+ ! adds absorbing boundary contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ endif
endif
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
@@ -663,9 +670,17 @@
potential_dot_dot_acoustic(iglob) &
- potential_dot_acoustic(iglob)*weight/cpl/rhol
elseif(SIMULATION_TYPE == 2) then
- potential_dot_dot_acoustic(iglob) = &
- potential_dot_dot_acoustic(iglob) &
- - b_absorb_acoustic_right(j,ib_right(ispecabs),NSTEP-it+1)
+ if(IS_BACKWARD_FIELD) then
+ ! adds (previously) stored contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - b_absorb_acoustic_right(j,ib_right(ispecabs),NSTEP-it+1)
+ else
+ ! adds absorbing boundary contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ endif
endif
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
@@ -702,9 +717,17 @@
potential_dot_dot_acoustic(iglob) &
- potential_dot_acoustic(iglob)*weight/cpl/rhol
elseif(SIMULATION_TYPE == 2) then
- potential_dot_dot_acoustic(iglob) = &
- potential_dot_dot_acoustic(iglob) &
- - b_absorb_acoustic_bottom(i,ib_bottom(ispecabs),NSTEP-it+1)
+ if(IS_BACKWARD_FIELD) then
+ ! adds (previously) stored contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - b_absorb_acoustic_bottom(i,ib_bottom(ispecabs),NSTEP-it+1)
+ else
+ ! adds absorbing boundary contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ endif
endif
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
@@ -741,9 +764,17 @@
potential_dot_dot_acoustic(iglob) &
- potential_dot_acoustic(iglob)*weight/cpl/rhol
elseif(SIMULATION_TYPE == 2) then
- potential_dot_dot_acoustic(iglob) = &
- potential_dot_dot_acoustic(iglob) &
- - b_absorb_acoustic_top(i,ib_top(ispecabs),NSTEP-it+1)
+ if(IS_BACKWARD_FIELD) then
+ ! adds (previously) stored contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - b_absorb_acoustic_top(i,ib_top(ispecabs),NSTEP-it+1)
+ else
+ ! adds absorbing boundary contribution
+ potential_dot_dot_acoustic(iglob) = &
+ potential_dot_dot_acoustic(iglob) &
+ - potential_dot_acoustic(iglob)*weight/cpl/rhol
+ endif
endif
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90 2011-09-07 17:01:35 UTC (rev 18877)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90 2011-09-07 17:23:20 UTC (rev 18878)
@@ -489,7 +489,7 @@
if(add_Bielak_conditions .and. initialfield) then
if (.not.over_critical_angle) then
call compute_Bielak_conditions(coord,iglob,nglob,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce(1), angleforce_refl, &
c_inc, c_refl, time_offset,f0)
traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
traction_z_t0 = mul_relaxed*(dxUz + dzUx)
@@ -580,7 +580,7 @@
if(add_Bielak_conditions .and. initialfield) then
if (.not.over_critical_angle) then
call compute_Bielak_conditions(coord,iglob,nglob,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce(1), angleforce_refl, &
c_inc, c_refl, time_offset,f0)
traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
traction_z_t0 = mul_relaxed*(dxUz + dzUx)
@@ -677,7 +677,7 @@
if(add_Bielak_conditions .and. initialfield) then
if (.not.over_critical_angle) then
call compute_Bielak_conditions(coord,iglob,nglob,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce(1), angleforce_refl, &
c_inc, c_refl, time_offset,f0)
traction_x_t0 = mul_relaxed*(dxUz + dzUx)
traction_z_t0 = lambdal_relaxed*dxUx + (lambdal_relaxed+2*mul_relaxed)*dzUz
@@ -773,7 +773,7 @@
! top or bottom edge, vertical normal vector
if(add_Bielak_conditions .and. initialfield) then
call compute_Bielak_conditions(coord,iglob,nglob,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce(1), angleforce_refl, &
c_inc, c_refl, time_offset,f0)
traction_x_t0 = mul_relaxed*(dxUz + dzUx)
traction_z_t0 = lambdal_relaxed*dxUx + (lambdal_relaxed+2*mul_relaxed)*dzUz
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-09-07 17:01:35 UTC (rev 18877)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-09-07 17:23:20 UTC (rev 18878)
@@ -195,7 +195,8 @@
!local parameters
integer :: it,i,j,iglob
real(kind=CUSTOM_REAL) :: t
- real(kind=CUSTOM_REAL) :: xigll, zigll
+ real(kind=CUSTOM_REAL), dimension(NGLLX) :: xigll
+ real(kind=CUSTOM_REAL), dimension(NGLLZ) :: zigll
real(kind=CUSTOM_REAL), dimension(NGLLX) :: hxi, hpxi
real(kind=CUSTOM_REAL), dimension(NGLLZ) :: hgamma, hpgamma
real(kind=CUSTOM_REAL) :: factor_noise, f0, aval, t0
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/paco_beyond_critical.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/paco_beyond_critical.f90 2011-09-07 17:01:35 UTC (rev 18877)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/paco_beyond_critical.f90 2011-09-07 17:23:20 UTC (rev 18878)
@@ -36,9 +36,9 @@
double precision, dimension(nbot,NSTEP_global) :: v0x_bot,v0z_bot, t0x_bot,t0z_bot
double precision, dimension(2,nglob) :: coord
- double precision, dimension(2,nglob) :: displ_elastic
- double precision, dimension(2,nglob) :: veloc_elastic
- double precision, dimension(2,nglob) :: accel_elastic
+ real(kind=CUSTOM_REAL), dimension(2,nglob) :: displ_elastic
+ real(kind=CUSTOM_REAL), dimension(2,nglob) :: veloc_elastic
+ real(kind=CUSTOM_REAL), dimension(2,nglob) :: accel_elastic
integer, dimension(:),allocatable :: local_pt
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90 2011-09-07 17:01:35 UTC (rev 18877)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/plotpost.F90 2011-09-07 17:23:20 UTC (rev 18878)
@@ -112,7 +112,8 @@
double precision coorg(NDIM,npgeo)
double precision, dimension(nrec) :: st_xval,st_zval
- integer numabs(nelemabs),codeabs(4,nelemabs)
+ integer numabs(nelemabs)
+ logical codeabs(4,nelemabs)
logical anyabs,coupled_acoustic_elastic,coupled_acoustic_poro,coupled_elastic_poro, &
any_acoustic,any_poroelastic,plot_lowerleft_corner_only
@@ -2205,7 +2206,7 @@
do iedge = 1,4
- if(codeabs(iedge,inum) /= 0) then
+ if(codeabs(iedge,inum)) then ! codeabs(:,:) is defined as "logical" in MAIN program
if(iedge == ITOP) then
ideb = 3
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-09-07 17:01:35 UTC (rev 18877)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-09-07 17:23:20 UTC (rev 18878)
@@ -1,4 +1,4 @@
-
+ program specfem2D
!========================================================================
!
! S P E C F E M 2 D Version 6 . 2
@@ -323,7 +323,6 @@
! not work because it would be discontinuous at such an interface and would
! therefore not be consistent with the basis functions.
- program specfem2D
implicit none
@@ -501,7 +500,7 @@
double precision, dimension(3):: bl_relaxed
double precision :: permlxx,permlxz,permlzz,invpermlxx,invpermlxz,invpermlzz,detk
! adjoint
- double precision, dimension(:), allocatable :: b_viscodampx,b_viscodampz
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: b_viscodampx,b_viscodampz
integer reclen
! for fluid/solid coupling and edge detection
@@ -791,6 +790,16 @@
integer :: ispecperio, ispecperio2, ispec2, i2, j2
integer :: iglob_target_to_replace, ispec3, i3, j3
+!<SU_FORMAT
+ integer(kind=4) :: header4(1)
+ integer(kind=2) :: header2(2)
+ equivalence(header2,header4)
+ integer(kind=4) :: r4head(60)
+ logical,parameter :: SU_FORMAT=.true.
+ character(len=512) :: filename
+ real(kind=4),dimension(:,:),allocatable :: adj_src_s
+!>SU_FORMAT
+
!<NOISE_TOMOGRAPHY
! NOISE_TOMOGRAPHY = 0 - turn noise tomography subroutines off, resulting in
! an earthquake simulation rather than a noise simulation
@@ -1769,34 +1778,68 @@
! compute source array for adjoint source
if(SIMULATION_TYPE == 2) then ! adjoint calculation
- nadj_rec_local = 0
- do irec = 1,nrec
- if(myrank == which_proc_receiver(irec))then
-! check that the source proc number is okay
- if(which_proc_receiver(irec) < 0 .or. which_proc_receiver(irec) > NPROC-1) &
- call exit_MPI('something is wrong with the source proc number in adjoint simulation')
- nadj_rec_local = nadj_rec_local + 1
- endif
- enddo
- if(ipass == 1) allocate(adj_sourcearray(NSTEP,3,NGLLX,NGLLZ))
- if (nadj_rec_local > 0 .and. ipass == 1) then
- allocate(adj_sourcearrays(nadj_rec_local,NSTEP,3,NGLLX,NGLLZ))
- else if (ipass == 1) then
- allocate(adj_sourcearrays(1,1,1,1,1))
+ if (.not. SU_FORMAT) then
+ nadj_rec_local = 0
+ do irec = 1,nrec
+ if(myrank == which_proc_receiver(irec))then
+ ! check that the source proc number is okay
+ if(which_proc_receiver(irec) < 0 .or. which_proc_receiver(irec) > NPROC-1) &
+ call exit_MPI('something is wrong with the source proc number in adjoint simulation')
+ nadj_rec_local = nadj_rec_local + 1
+ endif
+ enddo
+ if(ipass == 1) allocate(adj_sourcearray(NSTEP,3,NGLLX,NGLLZ))
+ if (nadj_rec_local > 0 .and. ipass == 1) then
+ allocate(adj_sourcearrays(nadj_rec_local,NSTEP,3,NGLLX,NGLLZ))
+ else if (ipass == 1) then
+ allocate(adj_sourcearrays(1,1,1,1,1))
+ endif
+
+ irec_local = 0
+ do irec = 1, nrec
+ ! compute only adjoint source arrays in the local proc
+ if(myrank == which_proc_receiver(irec))then
+ irec_local = irec_local + 1
+ adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
+ call compute_arrays_adj_source(adj_source_file, &
+ xi_receiver(irec), gamma_receiver(irec), &
+ adj_sourcearray, xigll,zigll,NSTEP)
+ adj_sourcearrays(irec_local,:,:,:,:) = adj_sourcearray(:,:,:,:)
+ endif
+ enddo
+ else
+ irec_local = 0
+ write(filename, "('./SEM/Ux_file_single.bin.adj')")
+ open(111,file=trim(filename),access='direct',recl=240+4*NSTEP,iostat = ios)
+ if (ios /= 0) call exit_MPI(' file '//trim(filename)//'does not exist')
+ write(filename, "('./SEM/Uz_file_single.bin.adj')")
+ open(113,file=trim(filename),access='direct',recl=240+4*NSTEP,iostat = ios)
+ if (ios /= 0) call exit_MPI(' file '//trim(filename)//'does not exist')
+
+ allocate(adj_src_s(NSTEP,3))
+
+ do irec = 1, nrec
+ irec_local = irec_local + 1
+ adj_sourcearray(:,:,:,:) = 0.0
+ read(111,rec=irec,iostat=ios) r4head, adj_src_s(:,1)
+ if (ios /= 0) call exit_MPI(' file '//trim(filename)//' read error')
+ read(113,rec=irec,iostat=ios) r4head, adj_src_s(:,3)
+ if (ios /= 0) call exit_MPI(' file '//trim(filename)//' read error (even if in SH case, this file is needed, although not used. so just make a duplication)')
+ header4=r4head(29)
+ !if (irec==1) print*, r4head(1),r4head(19),r4head(20),r4head(21),r4head(22),header2(2)
+ call lagrange_any(xi_receiver(irec),NGLLX,xigll,hxir,hpxir)
+ call lagrange_any(gamma_receiver(irec),NGLLZ,zigll,hgammar,hpgammar)
+ do k = 1, NGLLZ
+ do i = 1, NGLLX
+ adj_sourcearray(:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:)
+ enddo
+ enddo
+ adj_sourcearrays(irec_local,:,:,:,:) = adj_sourcearray(:,:,:,:)
+ enddo
+ close(111)
+ close(113)
+ deallocate(adj_src_s)
endif
-
- irec_local = 0
- do irec = 1, nrec
-! compute only adjoint source arrays in the local proc
- if(myrank == which_proc_receiver(irec))then
- irec_local = irec_local + 1
- adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
- call compute_arrays_adj_source(adj_source_file, &
- xi_receiver(irec), gamma_receiver(irec), &
- adj_sourcearray, xigll,zigll,NSTEP)
- adj_sourcearrays(irec_local,:,:,:,:) = adj_sourcearray(:,:,:,:)
- endif
- enddo
else if (ipass == 1) then
allocate(adj_sourcearrays(1,1,1,1,1))
endif
@@ -2906,7 +2949,7 @@
! call Paco's routine to compute in frequency and convert to time by Fourier transform
call paco_beyond_critical(coord,nglob,deltat,NSTEP,angleforce(1),&
- f0(1),cploc,csloc,TURN_ATTENUATION_ON,QKappa_attenuation,source_type(1),v0x_left,v0z_left,&
+ f0(1),cploc,csloc,TURN_ATTENUATION_ON,QKappa_attenuation(1),source_type(1),v0x_left,v0z_left,&
v0x_right,v0z_right,v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,&
t0x_bot,t0z_bot,left_bound(1:count_left),right_bound(1:count_right),bot_bound(1:count_bottom)&
,count_left,count_right,count_bottom,displ_paco,veloc_paco,accel_paco)
@@ -3994,7 +4037,7 @@
SIMULATION_TYPE,SAVE_FORWARD,nspec_left,nspec_right,&
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
b_absorb_acoustic_left,b_absorb_acoustic_right, &
- b_absorb_acoustic_bottom,b_absorb_acoustic_top)
+ b_absorb_acoustic_bottom,b_absorb_acoustic_top,.false.)
if( SIMULATION_TYPE == 2 ) then
call compute_forces_acoustic_2(nglob,nspec,nelemabs,numat,it,NSTEP, &
anyabs,assign_external_model,ibool,kmato,numabs, &
@@ -4008,7 +4051,7 @@
SIMULATION_TYPE,SAVE_FORWARD,nspec_left,nspec_right,&
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
b_absorb_acoustic_left,b_absorb_acoustic_right, &
- b_absorb_acoustic_bottom,b_absorb_acoustic_top)
+ b_absorb_acoustic_bottom,b_absorb_acoustic_top,.true.)
endif
@@ -4080,6 +4123,11 @@
if(SIMULATION_TYPE == 2) then
b_displ_x = b_displ_elastic(1,iglob)
b_displ_z = b_displ_elastic(3,iglob)
+ !<YANGL
+ ! new definition of adjoint displacement and adjoint potential
+ displ_x = accel_elastic(1,iglob)
+ displ_z = accel_elastic(3,iglob)
+ !>YANGL
endif
! get point values for the acoustic side
@@ -4301,7 +4349,7 @@
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_rec(irec))
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
+ + adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
enddo
enddo
endif ! if element acoustic
@@ -4512,6 +4560,10 @@
pressure = - potential_dot_dot_acoustic(iglob)
if(SIMULATION_TYPE == 2) then
b_pressure = - b_potential_dot_dot_acoustic(iglob)
+ !<YANGL
+ ! new definition of adjoint displacement and adjoint potential
+ pressure = potential_acoustic(iglob)
+ !>YANGL
endif
! get point values for the elastic side
ii2 = ivalue(ipoin1D,iedge_elastic)
@@ -6031,8 +6083,13 @@
do j = 1, NGLLZ
do i = 1, NGLLX
iglob = ibool(i,j,ispec)
- kappal_ac_global(iglob) = poroelastcoef(3,1,kmato(ispec))
- rhol_ac_global(iglob) = density(1,kmato(ispec))
+ if (.not. assign_external_model) then
+ kappal_ac_global(iglob) = poroelastcoef(3,1,kmato(ispec))
+ rhol_ac_global(iglob) = density(1,kmato(ispec))
+ else
+ rhol_ac_global(iglob) = rhoext(i,j,ispec)
+ kappal_ac_global(iglob) = rhoext(i,j,ispec)*vpext(i,j,ispec)*vpext(i,j,ispec)
+ endif
! calcul the displacement by computing the gradient of potential / rho
! and calcul the acceleration by computing the gradient of potential_dot_dot / rho
@@ -6042,11 +6099,13 @@
b_tempx2l = ZERO
do k = 1,NGLLX
! derivative along x
- tempx1l = tempx1l + potential_dot_dot_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
+ !tempx1l = tempx1l + potential_dot_dot_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
+ tempx1l = tempx1l + potential_acoustic(ibool(k,j,ispec))*hprime_xx(i,k) !!! YANGL
b_tempx1l = b_tempx1l + b_potential_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
bb_tempx1l = bb_tempx1l + b_potential_dot_dot_acoustic(ibool(k,j,ispec))*hprime_xx(i,k)
! derivative along z
- tempx2l = tempx2l + potential_dot_dot_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
+ !tempx2l = tempx2l + potential_dot_dot_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
+ tempx2l = tempx2l + potential_acoustic(ibool(i,k,ispec))*hprime_zz(j,k) !!! YANGL
b_tempx2l = b_tempx2l + b_potential_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
bb_tempx2l = bb_tempx2l + b_potential_dot_dot_acoustic(ibool(i,k,ispec))*hprime_zz(j,k)
enddo
@@ -6076,12 +6135,21 @@
do j = 1, NGLLZ
do i = 1, NGLLX
iglob = ibool(i,j,ispec)
- rho_ac_kl(i,j,ispec) = rho_ac_kl(i,j,ispec) - rhol_ac_global(iglob) * &
+ !<YANGL
+ !!!! old expression (from elastic kernels)
+ !!!rho_ac_kl(i,j,ispec) = rho_ac_kl(i,j,ispec) - rhol_ac_global(iglob) * &
+ !!! dot_product(accel_ac(:,iglob),b_displ_ac(:,iglob)) * deltat
+ !!!kappa_ac_kl(i,j,ispec) = kappa_ac_kl(i,j,ispec) - kappal_ac_global(iglob) * &
+ !!! potential_dot_dot_acoustic(iglob)/kappal_ac_global(iglob) * &
+ !!! b_potential_dot_dot_acoustic(iglob)/kappal_ac_global(iglob)&
+ !!! * deltat
+ !!!! new expression (from PDE-constrained optimization, coupling terms changed as well)
+ rho_ac_kl(i,j,ispec) = rho_ac_kl(i,j,ispec) + rhol_ac_global(iglob) * &
dot_product(accel_ac(:,iglob),b_displ_ac(:,iglob)) * deltat
- kappa_ac_kl(i,j,ispec) = kappa_ac_kl(i,j,ispec) - kappal_ac_global(iglob) * &
- potential_dot_dot_acoustic(iglob)/kappal_ac_global(iglob) * &
- b_potential_dot_dot_acoustic(iglob)/kappal_ac_global(iglob)&
- * deltat
+ kappa_ac_kl(i,j,ispec) = kappa_ac_kl(i,j,ispec) + kappal_ac_global(iglob) * &
+ potential_acoustic(iglob)/kappal_ac_global(iglob) * &
+ b_potential_dot_dot_acoustic(iglob)/kappal_ac_global(iglob) * deltat
+ !>YANGL
!
rhop_ac_kl(i,j,ispec) = rho_ac_kl(i,j,ispec) + kappa_ac_kl(i,j,ispec)
alpha_ac_kl(i,j,ispec) = TWO * kappa_ac_kl(i,j,ispec)
@@ -6103,10 +6171,17 @@
do j = 1, NGLLZ
do i = 1, NGLLX
iglob = ibool(i,j,ispec)
- mul_global(iglob) = poroelastcoef(2,1,kmato(ispec))
- kappal_global(iglob) = poroelastcoef(3,1,kmato(ispec)) &
- - 4._CUSTOM_REAL*mul_global(iglob)/3._CUSTOM_REAL
- rhol_global(iglob) = density(1,kmato(ispec))
+ if (.not. assign_external_model) then
+ mul_global(iglob) = poroelastcoef(2,1,kmato(ispec))
+ kappal_global(iglob) = poroelastcoef(3,1,kmato(ispec)) &
+ - 4._CUSTOM_REAL*mul_global(iglob)/3._CUSTOM_REAL
+ rhol_global(iglob) = density(1,kmato(ispec))
+ else
+ rhol_global(iglob) = rhoext(i,j,ispec)
+ mul_global(iglob) = rhoext(i,j,ispec)*vsext(i,j,ispec)*vsext(i,j,ispec)
+ kappal_global(iglob) = rhoext(i,j,ispec)*vpext(i,j,ispec)*vpext(i,j,ispec) &
+ -4._CUSTOM_REAL*mul_global(iglob)/3._CUSTOM_REAL
+ endif
rho_kl(i,j,ispec) = rho_kl(i,j,ispec) - rhol_global(iglob) * rho_k(iglob) * deltat
mu_kl(i,j,ispec) = mu_kl(i,j,ispec) - TWO * mul_global(iglob) * mu_k(iglob) * deltat
@@ -6339,9 +6414,10 @@
iglob = ibool(i,j,ispec)
xx = coord(1,iglob)
zz = coord(2,iglob)
- write(95,'(5e11.3)')xx,zz,rho_ac_kl(i,j,ispec),kappa_ac_kl(i,j,ispec)
- write(96,'(5e11.3)')rhorho_ac_hessian_final1(i,j,ispec), rhorho_ac_hessian_final2(i,j,ispec),&
- rhop_ac_kl(i,j,ispec),alpha_ac_kl(i,j,ispec)
+ write(95,'(4e15.5e4)')xx,zz,rho_ac_kl(i,j,ispec),kappa_ac_kl(i,j,ispec)
+ write(96,'(4e15.5e4)')xx,zz,rhop_ac_kl(i,j,ispec),alpha_ac_kl(i,j,ispec)
+ !write(96,'(4e15.5e4)')rhorho_ac_hessian_final1(i,j,ispec), rhorho_ac_hessian_final2(i,j,ispec),&
+ ! rhop_ac_kl(i,j,ispec),alpha_ac_kl(i,j,ispec)
enddo
enddo
enddo
@@ -6356,9 +6432,9 @@
iglob = ibool(i,j,ispec)
xx = coord(1,iglob)
zz = coord(2,iglob)
- write(97,'(5e11.3)')xx,zz,rho_kl(i,j,ispec),kappa_kl(i,j,ispec),mu_kl(i,j,ispec)
- write(98,'(5e11.3)')xx,zz,rhop_kl(i,j,ispec),alpha_kl(i,j,ispec),beta_kl(i,j,ispec)
- !write(98,'(5e11.3)')rhorho_el_hessian_final1(i,j,ispec), rhorho_el_hessian_final2(i,j,ispec),&
+ write(97,'(5e15.5e4)')xx,zz,rho_kl(i,j,ispec),kappa_kl(i,j,ispec),mu_kl(i,j,ispec)
+ write(98,'(5e15.5e4)')xx,zz,rhop_kl(i,j,ispec),alpha_kl(i,j,ispec),beta_kl(i,j,ispec)
+ !write(98,'(5e15.5e4)')rhorho_el_hessian_final1(i,j,ispec), rhorho_el_hessian_final2(i,j,ispec),&
! rhop_kl(i,j,ispec),alpha_kl(i,j,ispec),beta_kl(i,j,ispec)
enddo
enddo
@@ -6689,7 +6765,8 @@
if(.not. GENERATE_PARAVER_TRACES) &
call write_seismograms(sisux,sisuz,siscurl,station_name,network_name,NSTEP, &
nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
- NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current,p_sv)
+ NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current,p_sv,&
+ st_zval,x_source(1),z_source(1))
seismo_offset = seismo_offset + seismo_current
seismo_current = 0
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90 2011-09-07 17:01:35 UTC (rev 18877)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/write_seismograms.F90 2011-09-07 17:23:20 UTC (rev 18878)
@@ -47,7 +47,7 @@
subroutine write_seismograms(sisux,sisuz,siscurl,station_name,network_name, &
NSTEP,nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current,p_sv &
- )
+ ,st_zval,x_source,z_source)
implicit none
@@ -67,7 +67,7 @@
double precision, dimension(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc), intent(in) :: sisux,sisuz,siscurl
- double precision st_xval(nrec)
+ double precision :: st_xval(nrec)
character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name
character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name
@@ -87,6 +87,13 @@
integer :: irecloc
+!<SU_FORMAT
+ double precision :: st_zval(nrec),x_source,z_source
+ integer(kind=4) :: header4(1)
+ integer(kind=2) :: header2(2)
+ equivalence(header2,header4)
+!>SU_FORMAT
+
#ifdef USE_MPI
integer :: ierror
integer, dimension(MPI_STATUS_SIZE) :: status
@@ -151,9 +158,9 @@
open(unit=11,file='OUTPUT_FILES/Curl_file_double.bin',status='unknown')
close(11,status='delete')
- endif
+ endif
- if ( myrank == 0 ) then
+ if ( myrank == 0 ) then
! write the new files
if(seismotype == 4 .or. seismotype == 6) then
@@ -218,86 +225,123 @@
call MPI_RECV(buffer_binary(1,3),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(.not. SU_FORMAT) then
+ ! write trace
+ do iorientation = 1,number_of_components
- if(iorientation == 1) then
- chn = 'BXX'
- else if(iorientation == 2) then
- chn = 'BXZ'
- else if(iorientation == 3) then
- chn = 'cur'
- else
- call exit_MPI('incorrect channel value')
- endif
+ if(iorientation == 1) then
+ chn = 'BXX'
+ else if(iorientation == 2) then
+ chn = 'BXZ'
+ else if(iorientation == 3) then
+ chn = 'cur'
+ else
+ call exit_MPI('incorrect channel value')
+ endif
- ! in case of pressure, use different abbreviation
- if(seismotype == 4 .or. seismotype == 6) chn = 'PRE'
- ! in case of SH (membrane) waves, use different abbreviation
- if(.not.p_sv) chn = 'BXY'
+ ! in case of pressure, use different abbreviation
+ if(seismotype == 4 .or. seismotype == 6) chn = 'PRE'
+ ! in case of SH (membrane) waves, use different abbreviation
+ if(.not.p_sv) chn = 'BXY'
- ! 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))
- length_network_name = len_trim(network_name(irec))
+ ! 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))
+ length_network_name = len_trim(network_name(irec))
- ! 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
- call exit_MPI('wrong length of network name')
- end if
+ ! 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
+ 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
+ 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
- ! here would result in a loss of accuracy when one later convolves
- ! the results with the source time function
- if ( seismo_offset == 0 ) then
- open(unit=11,file=sisname(1:len_trim(sisname)),status='unknown')
- close(11,status='delete')
- endif
- open(unit=11,file=sisname(1:len_trim(sisname)),status='unknown',position='append')
+ ! 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
+ ! here would result in a loss of accuracy when one later convolves
+ ! the results with the source time function
+ if ( seismo_offset == 0 ) then
+ open(unit=11,file=sisname(1:len_trim(sisname)),status='unknown')
+ 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
- if(iorientation == 1) then
- write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', &
- sngl(buffer_binary(isample,iorientation))
- else
- write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', &
- sngl(buffer_binary(isample,iorientation))
- endif
- enddo
+ ! 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
+ if(iorientation == 1) then
+ write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', &
+ sngl(buffer_binary(isample,iorientation))
+ else
+ write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', &
+ sngl(buffer_binary(isample,iorientation))
+ endif
+ enddo
- close(11)
- end do
+ close(11)
+ end do
+ ! write binary seismogram
+ do isample = 1, seismo_current
+ write(12,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,1))
+ write(13,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,1)
+ if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then
+ write(14,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2))
+ write(15,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,2)
+ end if
+ if ( seismotype == 5 ) then
+ write(16,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,3))
+ write(17,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,3)
+ end if
+ enddo
+ else
+ if (seismo_offset==0) then
+ ! write SU headers (refer to Seismic Unix for details)
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+1) irec ! receiver ID
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+10) NINT(st_xval(irec)-x_source) ! offset
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+19) NINT(x_source) ! source location xs
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+20) NINT(z_source) ! source location zs
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+21) NINT(st_xval(irec)) ! receiver location xr
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+22) NINT(st_zval(irec)) ! receiver location zr
+ if (nrec>1) write(12,rec=(irec-1)*60+(irec-1)*NSTEP+48) SNGL(st_xval(2)-st_xval(1)) ! receiver interval
+ header2(1)=0 ! dummy
+ header2(2)=NSTEP
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+29) header4 ! equivalence(header4(1),header2(2))
+ header2(1)=NINT(deltat*1.0d6) ! deltat (unit: 10^{-6} second)
+ header2(2)=0 ! dummy
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+30) header4 ! equivalence(header4(1),header2(2))
+ if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then
+ ! headers
+ if (seismo_offset==0) then
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+1) irec
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+10) NINT(st_xval(irec)-x_source)
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+19) NINT(x_source)
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+20) NINT(z_source)
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+21) NINT(st_xval(irec))
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+22) NINT(st_zval(irec))
+ if(nrec>1) write(14,rec=(irec-1)*60+(irec-1)*NSTEP+48) SNGL(st_xval(2)-st_xval(1))
+ header2(1)=0 ! dummy
+ header2(2)=NSTEP
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+29) header4
+ header2(1)=NINT(deltat*1.0d6)
+ header2(2)=0 ! dummy
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+30) header4
+ end if
+ endif
+ endif
+ ! the "60" in the following corresponds to 240 bytes header (note the reclength is 4 bytes)
+ write(12,rec=irec*60+(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,1))
+ write(14,rec=irec*60+(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2))
+ endif
-! write binary seismogram
- do isample = 1, seismo_current
- write(12,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,1))
- write(13,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,1)
- if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then
- write(14,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2))
- write(15,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,2)
- end if
- if ( seismotype == 5 ) then
- write(16,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,3))
- write(17,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,3)
- end if
- enddo
#ifdef USE_MPI
-
else
if ( which_proc_receiver(irec) == myrank ) then
irecloc = irecloc + 1
@@ -309,7 +353,6 @@
call MPI_SEND(siscurl(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
end if
end if
-
#endif
end if
More information about the CIG-COMMITS
mailing list