[cig-commits] r8586 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:58:42 PST 2007
Author: walter
Date: 2007-12-07 15:58:41 -0800 (Fri, 07 Dec 2007)
New Revision: 8586
Modified:
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
Log:
data in the seismograms between two NTSTEP_BETWEEN_OUTPUT_INFO are now added to the binary seismograms, and appended to the ascii seismograms. The full seismograms are no longer written every NTSTEP_BETWEEN_OUTPUT_INFO; arrays for storing seismograms are now of size NTSTEP_BETWEEN_OUTPUT_INFO instead of NSTEP, and it runs faster when there are lots of receivers. This rev is slower when there are only a few receivers though.
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-09-13 13:34:46 UTC (rev 8585)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-07 23:58:41 UTC (rev 8586)
@@ -110,6 +110,7 @@
! for seismograms
double precision, dimension(:,:), allocatable :: sisux,sisuz
+ integer :: seismo_offset, seismo_current
! vector field in an element
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLX) :: vector_field_element
@@ -171,7 +172,7 @@
double precision :: vpmin,vpmax
- integer :: colors,numbers,subsamp,imagetype,NTSTEP_BETWEEN_OUTPUT_INFO,nrec,seismotype
+ integer :: colors,numbers,subsamp,imagetype,NTSTEP_BETWEEN_OUTPUT_INFO,NTSTEP_BETWEEN_OUTPUT_SEISMO,nrec,seismotype
integer :: numat,ngnod,nspec,pointsdisp,nelemabs,nelem_acoustic_surface,ispecabs
logical interpol,meshvect,modelvect,boundvect,assign_external_model,initialfield, &
@@ -374,6 +375,7 @@
read(IIN,"(a80)") datlin
read(IIN,*) NTSTEP_BETWEEN_OUTPUT_INFO
+ NTSTEP_BETWEEN_OUTPUT_SEISMO = NTSTEP_BETWEEN_OUTPUT_INFO
read(IIN,"(a80)") datlin
read(IIN,*) output_postscript_snapshot,output_color_image,colors,numbers
@@ -936,8 +938,8 @@
xi_receiver,gamma_receiver,station_name,network_name,x_source,z_source,coorg,knods,ngnod,npgeo)
! allocate seismogram arrays
- allocate(sisux(NSTEP,nrecloc))
- allocate(sisuz(NSTEP,nrecloc))
+ allocate(sisux(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
+ allocate(sisuz(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
! check if acoustic receiver is exactly on the free surface because pressure is zero there
do ispec_acoustic_surface = 1,nelem_acoustic_surface
@@ -1741,12 +1743,20 @@
#endif
+! initialize variables for writing seismograms
+ seismo_offset = 0
+ seismo_current = 0
+
+
! *********************************************************
! ************* MAIN LOOP OVER THE TIME STEPS *************
! *********************************************************
do it = 1,NSTEP
+! update position in seismograms
+ seismo_current = seismo_current + 1
+
! compute current time
time = (it-1)*deltat
@@ -2169,11 +2179,11 @@
! rotate seismogram components if needed, except if recording pressure, which is a scalar
if(seismotype /= 4) then
- sisux(it,irecloc) = cosrot*valux + sinrot*valuz
- sisuz(it,irecloc) = - sinrot*valux + cosrot*valuz
+ sisux(seismo_current,irecloc) = cosrot*valux + sinrot*valuz
+ sisuz(seismo_current,irecloc) = - sinrot*valux + cosrot*valuz
else
- sisux(it,irecloc) = valux
- sisuz(it,irecloc) = ZERO
+ sisux(seismo_current,irecloc) = valux
+ sisuz(seismo_current,irecloc) = ZERO
endif
enddo
@@ -2352,7 +2362,11 @@
!---- save temporary or final seismograms
call write_seismograms(sisux,sisuz,station_name,network_name,NSTEP, &
- nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,it,t0)
+ nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,it,t0, &
+ NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current &
+ )
+ seismo_offset = seismo_offset + seismo_current
+ seismo_current = 0
! count elapsed wall-clock time
call date_and_time(datein,timein,zone,time_values)
Modified: seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2007-09-13 13:34:46 UTC (rev 8585)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2007-12-07 23:58:41 UTC (rev 8586)
@@ -14,7 +14,9 @@
! write seismograms to text files
subroutine write_seismograms(sisux,sisuz,station_name,network_name, &
- NSTEP,nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,it,t0)
+ NSTEP,nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,it,t0, &
+ NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current &
+ )
implicit none
@@ -24,12 +26,14 @@
#endif
integer :: nrec,NSTEP,it,seismotype
+ 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
- double precision, dimension(NSTEP,nrecloc), intent(in) :: sisux,sisuz
+ double precision, dimension(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc), intent(in) :: sisux,sisuz
double precision st_xval(nrec)
@@ -81,10 +85,10 @@
number_of_components = NDIM
endif
- allocate(buffer_binary(NSTEP,number_of_components))
+ allocate(buffer_binary(NTSTEP_BETWEEN_OUTPUT_SEISMO,number_of_components))
- if ( myrank == 0 ) then
+ if ( myrank == 0 .and. seismo_offset == 0 ) then
! delete the old files
open(unit=11,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown')
@@ -105,23 +109,27 @@
open(unit=11,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown')
close(11,status='delete')
+ endif
+
+ if ( myrank == 0 ) then
+
! write the new files
if(seismotype == 4) then
- open(unit=12,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown',access='direct',recl=4*NSTEP)
+ open(unit=12,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown',access='direct',recl=4)
else
- open(unit=12,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown',access='direct',recl=4*NSTEP)
+ open(unit=12,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown',access='direct',recl=4)
endif
if(seismotype == 4) then
- open(unit=13,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown',access='direct',recl=8*NSTEP)
+ open(unit=13,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown',access='direct',recl=8)
else
- open(unit=13,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown',access='direct',recl=8*NSTEP)
+ open(unit=13,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown',access='direct',recl=8)
endif
! no Z component seismogram if pressurs
if(seismotype /= 4) then
- open(unit=14,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown',access='direct',recl=4*NSTEP)
- open(unit=15,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown',access='direct',recl=8*NSTEP)
+ 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
@@ -142,10 +150,10 @@
#ifdef USE_MPI
else
- call MPI_RECV(buffer_binary(1,1),NSTEP,MPI_DOUBLE_PRECISION,&
+ call MPI_RECV(buffer_binary(1,1),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,&
which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
if ( number_of_components == 2 ) then
- call MPI_RECV(buffer_binary(1,2),NSTEP,MPI_DOUBLE_PRECISION,&
+ 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
@@ -188,15 +196,19 @@
! 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
- open(unit=11,file=sisname(1:len_trim(sisname)),status='unknown')
+ 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,min(it,NSTEP)
+ do isample = 1,seismo_current
if(iorientation == 1) then
- write(11,*) sngl(dble(isample-1)*deltat - t0),' ',sngl(buffer_binary(isample,iorientation))
+ write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ',sngl(buffer_binary(isample,iorientation))
else
- write(11,*) sngl(dble(isample-1)*deltat - t0),' ',sngl(buffer_binary(isample,iorientation))
+ write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ',sngl(buffer_binary(isample,iorientation))
endif
enddo
@@ -204,21 +216,22 @@
end do
! write binary seismogram
- write(12,rec=irec) sngl(buffer_binary(:,1))
- write(13,rec=irec) buffer_binary(:,1)
+ 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 ) then
- write(14,rec=irec) sngl(buffer_binary(:,2))
- write(15,rec=irec) buffer_binary(:,2)
+ 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
-
+ enddo
#ifdef USE_MPI
else
if ( which_proc_receiver(irec) == myrank ) then
irecloc = irecloc + 1
- call MPI_SEND(sisux(1,irecloc),NSTEP,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
+ 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),NSTEP,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
+ call MPI_SEND(sisuz(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
end if
end if
More information about the cig-commits
mailing list