[cig-commits] r12653 - seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat Aug 16 08:29:47 PDT 2008
Author: dkomati1
Date: 2008-08-16 08:29:47 -0700 (Sat, 16 Aug 2008)
New Revision: 12653
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_seismograms.F90
Log:
moved rotation of seismograms from X,Y,Z to Earth,North,Vertical to write_seismograms subroutine in order to improve performance
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90 2008-08-16 14:03:38 UTC (rev 12652)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90 2008-08-16 15:29:47 UTC (rev 12653)
@@ -342,8 +342,8 @@
integer :: seismo_offset, seismo_current
#ifdef USE_MPI
integer :: nit_written
- double precision :: uxd, uyd, uzd
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: seismograms
+ double precision :: uxd,uyd,uzd
+ double precision, dimension(:,:), allocatable :: uxdstore,uydstore,uzdstore
#endif
integer :: i,j,k,ispec,iglob,iglob_mantle,iglob_inner_core
@@ -811,7 +811,7 @@
call exit_MPI(myrank, 'SIMULATION_TYPE could be only 1, 2, or 3')
if (SIMULATION_TYPE /= 1 .and. NSOURCES > 999999) &
- call exit_MPI(myrank, 'for adjoint simulations, NSOURCES <= 999999, if you need more change i6.6 in write_seismograms.f90')
+ call exit_MPI(myrank, 'for adjoint simulations, NSOURCES <= 999999, if you need more change i6.6 in write_seismograms.F90')
if (ATTENUATION_VAL .or. SIMULATION_TYPE /= 1 .or. SAVE_FORWARD .or. (MOVIE_VOLUME .and. SIMULATION_TYPE /= 3)) then
COMPUTE_AND_STORE_STRAIN = .true.
@@ -1548,13 +1548,17 @@
#ifdef USE_MPI
! allocate seismogram array
if (nrec_local > 0) then
- allocate(seismograms(NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
- if (ier /= 0 ) then
- print *,"ABORTING can not allocate in specfem3D while allocating seismograms ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
+ allocate(uxdstore(nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
+ allocate(uydstore(nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
+ allocate(uzdstore(nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
+ if (ier /= 0) then
+ print *,"ABORTING can not allocate in specfem3D while allocating seismograms ier=",ier
+ call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
+ endif
! initialize seismograms
- seismograms(:,:,:) = 0._CUSTOM_REAL
+ uxdstore(:,:) = 0.d0
+ uydstore(:,:) = 0.d0
+ uzdstore(:,:) = 0.d0
nit_written = 0
endif
#endif
@@ -2473,15 +2477,11 @@
endif ! of if FASTER_RECEIVERS_POINTS_ONLY
-! store North, East and Vertical components
-! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- seismograms(:,irec_local,seismo_current) = sngl(scale_displ*(nu(:,1,irec)*uxd + &
- nu(:,2,irec)*uyd + nu(:,3,irec)*uzd))
- else
- seismograms(:,irec_local,seismo_current) = scale_displ*(nu(:,1,irec)*uxd + &
- nu(:,2,irec)*uyd + nu(:,3,irec)*uzd)
- endif
+! store X, Y and Z components for now
+! will be converted to North, East and Vertical components later
+ uxdstore(irec_local,seismo_current) = uxd
+ uydstore(irec_local,seismo_current) = uyd
+ uzdstore(irec_local,seismo_current) = uzd
enddo
@@ -2490,7 +2490,7 @@
! write the current or final seismograms
if(seismo_current == NTSTEP_BETWEEN_OUTPUT_SEISMOS .or. it == it_end) then
- call write_seismograms(myrank,seismograms,number_receiver_global,station_name, &
+ call write_seismograms(myrank,uxdstore,uydstore,uzdstore,number_receiver_global,station_name, &
network_name,stlat,stlon,stele,nrec,nrec_local,DT,t0,it_end, &
yr_SAC,jda_SAC,ho_SAC,mi_SAC,sec_SAC,t_cmt_SAC, &
elat_SAC,elon_SAC,depth_SAC,mb_SAC,ename_SAC,cmt_lat_SAC,cmt_lon_SAC,&
@@ -2498,7 +2498,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM, &
OUTPUT_SEISMOS_SAC_BINARY,ROTATE_SEISMOGRAMS_RT,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
seismo_offset,seismo_current,WRITE_SEISMOGRAMS_BY_MASTER, &
- SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,one_seismogram)
+ SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,one_seismogram,scale_displ,nu)
if(myrank==0) then
write(IMAIN,*)
@@ -2506,12 +2506,14 @@
write(IMAIN,*)
endif
-! prepare to shift to the next interval to store seismograms
+! prepare to shift to the next time interval to store seismograms
seismo_offset = seismo_offset + seismo_current
seismo_current = 0
-! clean seismogram array
- seismograms(:,:,:) = 0._CUSTOM_REAL
+! clean seismograms for the next time interval
+ uxdstore(:,:) = 0.d0
+ uydstore(:,:) = 0.d0
+ uzdstore(:,:) = 0.d0
endif
#endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_seismograms.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_seismograms.F90 2008-08-16 14:03:38 UTC (rev 12652)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_seismograms.F90 2008-08-16 15:29:47 UTC (rev 12653)
@@ -26,7 +26,7 @@
!=====================================================================
! write seismograms to files
- subroutine write_seismograms(myrank,seismograms,number_receiver_global,station_name, &
+ subroutine write_seismograms(myrank,uxdstore,uydstore,uzdstore,number_receiver_global,station_name, &
network_name,stlat,stlon,stele,nrec,nrec_local,DT,hdur,it_end, &
yr,jda,ho,mi,sec,t_cmt, &
elat,elon,depth,mb,ename,cmt_lat,cmt_lon, &
@@ -34,7 +34,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM, &
OUTPUT_SEISMOS_SAC_BINARY,ROTATE_SEISMOGRAMS_RT,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
seismo_offset,seismo_current,WRITE_SEISMOGRAMS_BY_MASTER,&
- SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,one_seismogram)
+ SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,one_seismogram,scale_displ,nu)
implicit none
@@ -55,8 +55,10 @@
integer :: seismo_offset, seismo_current, NTSTEP_BETWEEN_OUTPUT_SEISMOS
integer, dimension(nrec_local) :: number_receiver_global
+ double precision, dimension(nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS) :: uxdstore,uydstore,uzdstore
real(kind=CUSTOM_REAL), dimension(NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS) :: seismograms
- double precision hdur,DT
+ double precision, dimension(NDIM,NDIM,nrec) :: nu
+ double precision hdur,DT,scale_displ
character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name
character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name
@@ -70,7 +72,7 @@
! variables
integer :: iproc,sender,irec_local,irec,receiver,nrec_local_received,nrec_tot_found
- integer :: total_seismos,total_seismos_local
+ integer :: total_seismos,total_seismos_local,it_seismo,idimension
double precision :: write_time_begin,write_time
#ifdef USE_MPI
@@ -114,6 +116,30 @@
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
+! convert seismograms from X, Y, Z components to North, East and Vertical components
+! distinguish between single and double precision for reals
+ do it_seismo = 1,NTSTEP_BETWEEN_OUTPUT_SEISMOS
+ do irec_local = 1,nrec_local
+
+! get global number of that receiver
+ irec = number_receiver_global(irec_local)
+
+ do idimension = 1,NDIM
+ if(CUSTOM_REAL == SIZE_REAL) then
+ seismograms(idimension,irec_local,it_seismo) = &
+ sngl(scale_displ*(nu(idimension,1,irec)*uxdstore(irec_local,it_seismo) + &
+ nu(idimension,2,irec)*uydstore(irec_local,it_seismo) + &
+ nu(idimension,3,irec)*uzdstore(irec_local,it_seismo)))
+ else
+ seismograms(idimension,irec_local,it_seismo) = &
+ scale_displ*(nu(idimension,1,irec)*uxdstore(irec_local,it_seismo) + &
+ nu(idimension,2,irec)*uydstore(irec_local,it_seismo) + &
+ nu(idimension,3,irec)*uzdstore(irec_local,it_seismo))
+ endif
+ enddo
+ enddo
+ enddo
+
! all the processes write their local seismograms themselves
if(.not. WRITE_SEISMOGRAMS_BY_MASTER) then
@@ -334,7 +360,7 @@
real(kind=CUSTOM_REAL), dimension(5,NTSTEP_BETWEEN_OUTPUT_SEISMOS) :: seismogram_tmp
- integer myrank
+ integer myrank
double precision hdur,DT
character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name
@@ -434,8 +460,8 @@
endif
!do iorientation = 1,NDIM
- !do iorientation = 1,5 ! BS BS changed from 3 (NEZ) to 5 (NEZRT) components
- do iorientation = ior_start,ior_end ! BS BS changed according to ROTATE_SEISMOGRAMS_RT
+ !do iorientation = 1,5 ! changed from 3 (NEZ) to 5 (NEZRT) components
+ do iorientation = ior_start,ior_end ! changed according to ROTATE_SEISMOGRAMS_RT
if(iorientation == 1) then
chn = 'LHN'
@@ -451,9 +477,9 @@
call exit_MPI(myrank,'incorrect channel value')
endif
- if (iorientation == 4 .or. iorientation == 5) then ! LMU BS BS
+ if (iorientation == 4 .or. iorientation == 5) then
- ! BS BS calculate backazimuth needed to rotate East and North
+ ! calculate backazimuth needed to rotate East and North
! components to Radial and Transverse components
if (backaz>180.) then
@@ -467,7 +493,7 @@
cphi=cos(phi*pi/180)
sphi=sin(phi*pi/180)
- ! BS BS do the rotation of the components and put result in
+ ! rotate of the components and put result in
! new variable seismogram_tmp
if (iorientation == 4) then ! radial component
do isample = 1,seismo_current
More information about the cig-commits
mailing list