[cig-commits] r18164 - in seismo/3D/SPECFEM3D_GLOBE/trunk: setup src/auxiliaries src/create_header_file src/shared src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Fri Apr 1 21:05:46 PDT 2011
Author: danielpeter
Date: 2011-04-01 21:05:46 -0700 (Fri, 01 Apr 2011)
New Revision: 18164
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_surface.f90
Log:
updates create_movie_AVS_DX.f90; updates saving source and receiver vtk files; updates allocation checks and fixes compiler warnings; adds movie surface outputs for displacements if MOVIE_VOLUME_TYPE is set to 5
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2011-04-02 04:05:46 UTC (rev 18164)
@@ -165,7 +165,7 @@
!!-----------------------------------------------------------
!!
-!! Gauss-Lobatto-Legrendre resolution
+!! Gauss-Lobatto-Legendre resolution
!!
!!-----------------------------------------------------------
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90 2011-04-02 04:05:46 UTC (rev 18164)
@@ -91,20 +91,29 @@
include "constants.h"
+!---------------------
+! threshold and normalization parameters
+
+! single parameter to turn off all thresholding
+ logical, parameter :: ALL_THRESHOLD_OFF = .true.
+
! threshold in percent of the maximum below which we cut the amplitude
real(kind=CUSTOM_REAL), parameter :: THRESHOLD = 1._CUSTOM_REAL / 100._CUSTOM_REAL
! flag to apply non linear scaling to normalized norm of displacement
logical, parameter :: NONLINEAR_SCALING = .false.
- logical, parameter :: FIX_SCALING = .false. ! uses fixed max_value to normalize instead of max of current wavefield
- real,parameter:: MAX_VALUE = 6.77e-4
+! uses fixed max_value to normalize instead of max of current wavefield
+ logical, parameter :: FIX_SCALING = .true.
+ real,parameter:: MAX_VALUE = 1.0
! coefficient of power law used for non linear scaling
real(kind=CUSTOM_REAL), parameter :: POWER_SCALING = 0.30_CUSTOM_REAL
! flag to cut amplitude below a certain threshold
- logical, parameter :: APPLY_THRESHOLD = .true.
+ logical, parameter :: APPLY_THRESHOLD = .false.
+!---------------------
+
integer i,j,it
integer it1,it2
integer nspectot_AVS_max
@@ -398,64 +407,70 @@
enddo
-! compute min and max of data value to normalize
- min_field_current = minval(field_display(:))
- max_field_current = maxval(field_display(:))
+!---------------------------------
+! apply normalization and thresholding, if desired
-! make sure range is always symmetric and center is in zero
-! this assumption works only for fields that can be negative
-! would not work for norm of vector for instance
-! (we would lose half of the color palette if no negative values)
- max_absol = max(abs(min_field_current),abs(max_field_current))
- min_field_current = - max_absol
- max_field_current = + max_absol
+ if (.not. ALL_THRESHOLD_OFF) then
-! print minimum and maximum amplitude in current snapshot
- print *
- print *,'minimum amplitude in current snapshot = ',min_field_current
- print *,'maximum amplitude in current snapshot = ',max_field_current
- if( FIX_SCALING ) then
- print *,' to be normalized by : ',MAX_VALUE
- if( max_field_current > MAX_VALUE ) stop 'increase MAX_VALUE'
- endif
- print *
+ ! compute min and max of data value to normalize
+ min_field_current = minval(field_display(:))
+ max_field_current = maxval(field_display(:))
+ ! make sure range is always symmetric and center is in zero
+ ! this assumption works only for fields that can be negative
+ ! would not work for norm of vector for instance
+ ! (we would lose half of the color palette if no negative values)
+ max_absol = max(abs(min_field_current),abs(max_field_current))
+ min_field_current = - max_absol
+ max_field_current = + max_absol
+ ! print minimum and maximum amplitude in current snapshot
+ print *
+ print *,'minimum amplitude in current snapshot = ',min_field_current
+ print *,'maximum amplitude in current snapshot = ',max_field_current
+ if( FIX_SCALING ) then
+ print *,' to be normalized by : ',MAX_VALUE
+ if( max_field_current > MAX_VALUE ) stop 'increase MAX_VALUE'
+ endif
+ print *
-! normalize field to [0:1]
- print *,'normalizing... '
- if( FIX_SCALING ) then
- field_display(:) = (field_display(:) + MAX_VALUE) / (2.0*MAX_VALUE)
- else
- field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
- endif
-! rescale to [-1,1]
- field_display(:) = 2.*field_display(:) - 1.
+ ! normalize field to [0:1]
+ print *,'normalizing... '
+ if( FIX_SCALING ) then
+ field_display(:) = (field_display(:) + MAX_VALUE) / (2.0*MAX_VALUE)
+ else
+ field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
+ endif
+ ! rescale to [-1,1]
+ field_display(:) = 2.*field_display(:) - 1.
-! apply threshold to normalized field
- if(APPLY_THRESHOLD) then
- print *,'thresholding... '
- where(abs(field_display(:)) <= THRESHOLD) field_display = 0.
- endif
+ ! apply threshold to normalized field
+ if(APPLY_THRESHOLD) then
+ print *,'thresholding... '
+ where(abs(field_display(:)) <= THRESHOLD) field_display = 0.
+ endif
-! apply non linear scaling to normalized field if needed
- if(NONLINEAR_SCALING) then
- print *,'nonlinear scaling... '
- where(field_display(:) >= 0.)
- field_display = field_display ** POWER_SCALING
- elsewhere
- field_display = - abs(field_display) ** POWER_SCALING
- endwhere
- endif
+ ! apply non linear scaling to normalized field if needed
+ if(NONLINEAR_SCALING) then
+ print *,'nonlinear scaling... '
+ where(field_display(:) >= 0.)
+ field_display = field_display ** POWER_SCALING
+ elsewhere
+ field_display = - abs(field_display) ** POWER_SCALING
+ endwhere
+ endif
- print *,'color scaling... '
-! map back to [0,1]
- field_display(:) = (field_display(:) + 1.) / 2.
+ print *,'color scaling... '
+ ! map back to [0,1]
+ field_display(:) = (field_display(:) + 1.) / 2.
-! map field to [0:255] for AVS color scale
- field_display(:) = 255. * field_display(:)
+ ! map field to [0:255] for AVS color scale
+ field_display(:) = 255. * field_display(:)
+ endif
+!---------------------------------
+
! copy coordinate arrays since the sorting routine does not preserve them
print *,'sorting... '
xp_save(:) = xp(:)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90 2011-04-02 04:05:46 UTC (rev 18164)
@@ -215,9 +215,10 @@
print *,'time-stepping of the solver will be: ',DT
print *
if(MOVIE_SURFACE .or. MOVIE_VOLUME) then
- print *,'MOVIE_VOLUME:',MOVIE_VOLUME
+ print *,'MOVIE_VOLUME :',MOVIE_VOLUME
print *,'MOVIE_SURFACE:',MOVIE_SURFACE
print *,'Saving movie frames every',NTSTEP_BETWEEN_FRAMES
+ print *
endif
print *,'on NEC SX, make sure "loopcnt=" parameter'
! use fused loops on NEC SX
@@ -229,10 +230,14 @@
print *
print *,'size of static arrays per slice = ',static_memory_size/1073741824.d0,' GB'
print *
- print *,' (should be below and typically equal to 80% or 90%'
- print *,' of the memory installed per core)'
- print *,' (if significantly more, the job will not run by lack of memory)'
- print *,' (if significantly less, you waste a significant amount of memory)'
+ ! note: using less memory becomes only an issue if the code scaling is bad.
+ ! most users will run simulations with an executable using far less than 80% RAM per core
+ ! since they prefer having a faster computational time (and use a higher number of cores).
+ !print *,' (should be below and typically equal to 80% or 90%'
+ !print *,' of the memory installed per core)'
+ print *,' (should be below to 80% or 90% of the memory installed per core)'
+ print *,' (if significantly more, the job will not run by lack of memory )'
+ !print *,' (if significantly less, you waste a significant amount of memory)'
print *
print *,'size of static arrays for all slices = ',static_memory_size*dble(NPROCTOT)/1073741824.d0,' GB'
print *,' = ',static_memory_size*dble(NPROCTOT)/1099511627776.d0,' TB'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.f90 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.f90 2011-04-02 04:05:46 UTC (rev 18164)
@@ -261,10 +261,14 @@
write(IOUT,*) '!'
write(IOUT,*) '! size of static arrays per slice = ',static_memory_size/1073741824.d0,' GB'
write(IOUT,*) '!'
- write(IOUT,*) '! (should be below and typically equal to 80% or 90%'
- write(IOUT,*) '! of the memory installed per core)'
- write(IOUT,*) '! (if significantly more, the job will not run by lack of memory)'
- write(IOUT,*) '! (if significantly less, you waste a significant amount of memory)'
+ ! note: using less memory becomes only an issue if the code scaling is bad.
+ ! most users will run simulations with an executable using far less than 80% RAM per core
+ ! since they prefer having a faster computational time (and use a higher number of cores).
+ !write(IOUT,*) '! (should be below and typically equal to 80% or 90%'
+ !write(IOUT,*) '! of the memory installed per core)'
+ write(IOUT,*) '! (should be below to 80% or 90% of the memory installed per core)'
+ write(IOUT,*) '! (if significantly more, the job will not run by lack of memory )'
+ !write(IOUT,*) '! (if significantly less, you waste a significant amount of memory)'
write(IOUT,*) '!'
write(IOUT,*) '! size of static arrays for all slices = ',static_memory_size*dble(NPROCTOT)/1073741824.d0,' GB'
write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1099511627776.d0,' TB'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in 2011-04-02 04:05:46 UTC (rev 18164)
@@ -313,6 +313,9 @@
$O/compute_add_sources.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_add_sources.f90
${FCCOMPILE_NO_CHECK} -c -o $O/compute_add_sources.o ${FCFLAGS_f90} $S/compute_add_sources.f90
+$O/compute_boundary_kernel.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_boundary_kernel.f90
+ ${FCCOMPILE_CHECK} -c -o $O/compute_boundary_kernel.o ${FCFLAGS_f90} $S/compute_boundary_kernel.f90
+
$O/compute_coupling.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_coupling.f90
${FCCOMPILE_NO_CHECK} -c -o $O/compute_coupling.o ${FCFLAGS_f90} $S/compute_coupling.f90
@@ -334,6 +337,9 @@
$O/compute_forces_inner_core_Dev.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_forces_inner_core_Dev.f90
${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_inner_core_Dev.o ${FCFLAGS_f90} $S/compute_forces_inner_core_Dev.f90
+$O/compute_kernels.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_kernels.f90
+ ${FCCOMPILE_NO_CHECK} -c -o $O/compute_kernels.o ${FCFLAGS_f90} $S/compute_kernels.f90
+
$O/compute_seismograms.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/compute_seismograms.f90
${FCCOMPILE_NO_CHECK} -c -o $O/compute_seismograms.o ${FCFLAGS_f90} $S/compute_seismograms.f90
@@ -370,12 +376,6 @@
$O/compute_arrays_source.o: ${SETUP}/constants.h $S/compute_arrays_source.f90
${FCCOMPILE_CHECK} -c -o $O/compute_arrays_source.o ${FCFLAGS_f90} $S/compute_arrays_source.f90
-$O/compute_boundary_kernel.o: ${SETUP}/constants.h $S/compute_boundary_kernel.f90
- ${FCCOMPILE_CHECK} -c -o $O/compute_boundary_kernel.o ${FCFLAGS_f90} $S/compute_boundary_kernel.f90
-
-$O/compute_kernels.o: ${SETUP}/constants.h $S/compute_kernels.f90
- ${FCCOMPILE_CHECK} -c -o $O/compute_kernels.o ${FCFLAGS_f90} $S/compute_kernels.f90
-
$O/convert_time.o: $S/convert_time.f90
${FCCOMPILE_CHECK} -c -o $O/convert_time.o ${FCFLAGS_f90} $S/convert_time.f90
@@ -409,9 +409,15 @@
$O/check_simulation_stability.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/check_simulation_stability.f90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/check_simulation_stability.o ${FCFLAGS_f90} $S/check_simulation_stability.f90
+$O/fix_non_blocking_flags.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/fix_non_blocking_flags.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/fix_non_blocking_flags.o ${FCFLAGS_f90} $S/fix_non_blocking_flags.f90
+
$O/initialize_simulation.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/initialize_simulation.f90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/initialize_simulation.o ${FCFLAGS_f90} $S/initialize_simulation.f90
+$O/noise_tomography.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/noise_tomography.f90
+ ${MPIFCCOMPILE_NO_CHECK} -c -o $O/noise_tomography.o ${FCFLAGS_f90} $S/noise_tomography.f90
+
$O/prepare_timerun.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/prepare_timerun.f90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/prepare_timerun.o ${FCFLAGS_f90} $S/prepare_timerun.f90
@@ -424,17 +430,12 @@
$O/specfem3D.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/specfem3D.f90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/specfem3D.o ${FCFLAGS_f90} $S/specfem3D.f90
-$O/fix_non_blocking_flags.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/fix_non_blocking_flags.f90
- ${MPIFCCOMPILE_CHECK} -c -o $O/fix_non_blocking_flags.o ${FCFLAGS_f90} $S/fix_non_blocking_flags.f90
-
$O/write_movie_surface.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/write_movie_surface.f90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/write_movie_surface.o ${FCFLAGS_f90} $S/write_movie_surface.f90
$O/write_movie_volume.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/write_movie_volume.f90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/write_movie_volume.o ${FCFLAGS_f90} $S/write_movie_volume.f90
-$O/noise_tomography.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/noise_tomography.f90
- ${MPIFCCOMPILE_NO_CHECK} -c -o $O/noise_tomography.o ${FCFLAGS_f90} $S/noise_tomography.f90
##
## specfem3D - non-dependent on values from mesher here
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90 2011-04-02 04:05:46 UTC (rev 18164)
@@ -175,27 +175,27 @@
endif
! allocate memory for arrays using number of stations
- allocate(epidist(nrec))
- allocate(ix_initial_guess(nrec))
- allocate(iy_initial_guess(nrec))
- allocate(iz_initial_guess(nrec))
- allocate(x_target(nrec))
- allocate(y_target(nrec))
- allocate(z_target(nrec))
- allocate(x_found(nrec))
- allocate(y_found(nrec))
- allocate(z_found(nrec))
- allocate(final_distance(nrec))
+ allocate(epidist(nrec), &
+ ix_initial_guess(nrec), &
+ iy_initial_guess(nrec), &
+ iz_initial_guess(nrec), &
+ x_target(nrec), &
+ y_target(nrec), &
+ z_target(nrec), &
+ x_found(nrec), &
+ y_found(nrec), &
+ z_found(nrec), &
+ final_distance(nrec), &
+ ispec_selected_rec_all(nrec,0:NPROCTOT-1), &
+ xi_receiver_all(nrec,0:NPROCTOT-1), &
+ eta_receiver_all(nrec,0:NPROCTOT-1), &
+ gamma_receiver_all(nrec,0:NPROCTOT-1), &
+ x_found_all(nrec,0:NPROCTOT-1), &
+ y_found_all(nrec,0:NPROCTOT-1), &
+ z_found_all(nrec,0:NPROCTOT-1), &
+ final_distance_all(nrec,0:NPROCTOT-1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating temporary receiver arrays')
- allocate(ispec_selected_rec_all(nrec,0:NPROCTOT-1))
- allocate(xi_receiver_all(nrec,0:NPROCTOT-1))
- allocate(eta_receiver_all(nrec,0:NPROCTOT-1))
- allocate(gamma_receiver_all(nrec,0:NPROCTOT-1))
- allocate(x_found_all(nrec,0:NPROCTOT-1))
- allocate(y_found_all(nrec,0:NPROCTOT-1))
- allocate(z_found_all(nrec,0:NPROCTOT-1))
- allocate(final_distance_all(nrec,0:NPROCTOT-1))
-
! read that STATIONS file on the master
if(myrank == 0) then
call get_value_string(STATIONS, 'solver.STATIONS', rec_filename)
@@ -319,7 +319,8 @@
y_target(irec) = r0*sin(theta)*sin(phi)
z_target(irec) = r0*cos(theta)
- if (myrank == 0) write(IOVTK,*) sngl(x_target(irec)), sngl(y_target(irec)), sngl(z_target(irec))
+ ! would write out desired target locations of receivers
+ !if (myrank == 0) write(IOVTK,*) sngl(x_target(irec)), sngl(y_target(irec)), sngl(z_target(irec))
! examine top of the elements only (receivers always at the surface)
! k = NGLLZ
@@ -541,21 +542,28 @@
ispec_selected_rec_all(:,:) = -1
call MPI_GATHER(ispec_selected_rec,nrec,MPI_INTEGER,ispec_selected_rec_all,nrec,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
- call MPI_GATHER(xi_receiver,nrec,MPI_DOUBLE_PRECISION,xi_receiver_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
- call MPI_GATHER(eta_receiver,nrec,MPI_DOUBLE_PRECISION,eta_receiver_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
- call MPI_GATHER(gamma_receiver,nrec,MPI_DOUBLE_PRECISION,gamma_receiver_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
- call MPI_GATHER(final_distance,nrec,MPI_DOUBLE_PRECISION,final_distance_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
- call MPI_GATHER(x_found,nrec,MPI_DOUBLE_PRECISION,x_found_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
- call MPI_GATHER(y_found,nrec,MPI_DOUBLE_PRECISION,y_found_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
- call MPI_GATHER(z_found,nrec,MPI_DOUBLE_PRECISION,z_found_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(xi_receiver,nrec,MPI_DOUBLE_PRECISION,xi_receiver_all,nrec, &
+ MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(eta_receiver,nrec,MPI_DOUBLE_PRECISION,eta_receiver_all,nrec, &
+ MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(gamma_receiver,nrec,MPI_DOUBLE_PRECISION,gamma_receiver_all,nrec, &
+ MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(final_distance,nrec,MPI_DOUBLE_PRECISION,final_distance_all,nrec, &
+ MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(x_found,nrec,MPI_DOUBLE_PRECISION,x_found_all,nrec, &
+ MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(y_found,nrec,MPI_DOUBLE_PRECISION,y_found_all,nrec, &
+ MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(z_found,nrec,MPI_DOUBLE_PRECISION,z_found_all,nrec, &
+ MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
! this is executed by main process only
if(myrank == 0) then
-! check that the gather operation went well
+ ! 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
+ ! MPI loop on all the results to determine the best slice
islice_selected_rec(:) = -1
do irec = 1,nrec
distmin = HUGEVAL
@@ -591,8 +599,8 @@
write(IMAIN,*) ' at xi,eta,gamma coordinates = ',xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec)
endif
-! add warning if estimate is poor
-! (usually means receiver outside the mesh given by the user)
+ ! add warning if estimate is poor
+ ! (usually means receiver outside the mesh given by the user)
if(final_distance(irec) > THRESHOLD_EXCLUDE_STATION) then
write(IMAIN,*) 'station # ',irec,' ',station_name(irec),network_name(irec)
write(IMAIN,*) '*****************************************************************'
@@ -604,6 +612,7 @@
write(IMAIN,*) '*****************************************************************'
else
nrec_found = nrec_found + 1
+
islice_selected_rec_found(nrec_found) = islice_selected_rec(irec)
ispec_selected_rec_found(nrec_found) = ispec_selected_rec(irec)
xi_receiver_found(nrec_found) = xi_receiver(irec)
@@ -617,19 +626,22 @@
stbur_found(nrec_found) = stbur(irec)
nu_found(:,:,nrec_found) = nu(:,:,irec)
epidist_found(nrec_found) = epidist(irec)
+
+ ! writes out actual receiver location to vtk file
+ write(IOVTK,*) sngl(x_found(irec)), sngl(y_found(irec)), sngl(z_found(irec))
endif
enddo
-! compute maximal distance for all the receivers
+ ! compute maximal distance for all the receivers
final_distance_max = maxval(final_distance(:))
-! display maximum error for all the receivers
+ ! display maximum error for all the receivers
write(IMAIN,*)
write(IMAIN,*) 'maximum error in location of all the receivers: ',sngl(final_distance_max),' km'
-! add warning if estimate is poor
-! (usually means receiver outside the mesh given by the user)
+ ! add warning if estimate is poor
+ ! (usually means receiver outside the mesh given by the user)
if(final_distance_max > THRESHOLD_EXCLUDE_STATION) then
write(IMAIN,*)
write(IMAIN,*) '************************************************************'
@@ -679,9 +691,7 @@
close(27)
endif
-
-
-! elapsed time since beginning of mesh generation
+ ! elapsed time since beginning of mesh generation
tCPU = MPI_WTIME() - time_start
write(IMAIN,*)
write(IMAIN,*) 'Elapsed time for receiver detection in seconds = ',tCPU
@@ -710,7 +720,7 @@
call MPI_BCAST(stbur,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
call MPI_BCAST(nu,nrec*3*3,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! deallocate arrays
+ ! deallocate arrays
deallocate(epidist)
deallocate(ix_initial_guess)
deallocate(iy_initial_guess)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90 2011-04-02 04:05:46 UTC (rev 18164)
@@ -226,30 +226,24 @@
NSOURCES_SUBSET_current_size = min(NSOURCES_SUBSET_MAX, NSOURCES - isources_already_done)
! allocate arrays specific to each subset
- allocate(final_distance_source_subset(NSOURCES_SUBSET_current_size))
+ allocate(final_distance_source_subset(NSOURCES_SUBSET_current_size), &
+ ispec_selected_source_subset(NSOURCES_SUBSET_current_size), &
+ ispec_selected_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1), &
+ xi_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1), &
+ eta_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1), &
+ gamma_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1), &
+ final_distance_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1), &
+ x_found_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1), &
+ y_found_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1), &
+ z_found_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1), &
+ xi_source_subset(NSOURCES_SUBSET_current_size), &
+ eta_source_subset(NSOURCES_SUBSET_current_size), &
+ gamma_source_subset(NSOURCES_SUBSET_current_size), &
+ x_found_source(NSOURCES_SUBSET_current_size), &
+ y_found_source(NSOURCES_SUBSET_current_size), &
+ z_found_source(NSOURCES_SUBSET_current_size),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating temporary source arrays')
- allocate(ispec_selected_source_subset(NSOURCES_SUBSET_current_size))
-
- allocate(ispec_selected_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1))
-
- allocate(xi_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1))
- allocate(eta_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1))
- allocate(gamma_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1))
-
- allocate(final_distance_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1))
-
- allocate(x_found_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1))
- allocate(y_found_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1))
- allocate(z_found_source_all(NSOURCES_SUBSET_current_size,0:NPROCTOT-1))
-
- allocate(xi_source_subset(NSOURCES_SUBSET_current_size))
- allocate(eta_source_subset(NSOURCES_SUBSET_current_size))
- allocate(gamma_source_subset(NSOURCES_SUBSET_current_size))
-
- allocate(x_found_source(NSOURCES_SUBSET_current_size))
- allocate(y_found_source(NSOURCES_SUBSET_current_size))
- allocate(z_found_source(NSOURCES_SUBSET_current_size))
-
! make sure we clean the subset array before the gather
ispec_selected_source_subset(:) = 0
@@ -355,7 +349,8 @@
y_target_source = r_target_source*dsin(theta)*dsin(phi)
z_target_source = r_target_source*dcos(theta)
- if(myrank == 0) write(IOVTK,*) sngl(x_target_source),sngl(y_target_source),sngl(z_target_source)
+ ! would only output desired target locations
+ !if(myrank == 0) write(IOVTK,*) sngl(x_target_source),sngl(y_target_source),sngl(z_target_source)
! set distance to huge initial value
distmin = HUGEVAL
@@ -595,205 +590,210 @@
! this is executed by main process only
if(myrank == 0) then
-! check that the gather operation went well
+ ! check that the gather operation went well
if(minval(ispec_selected_source_all) <= 0) call exit_MPI(myrank,'gather operation failed for source')
-! loop on all the sources within subsets
+ ! loop on all the sources within subsets
do isource_in_this_subset = 1,NSOURCES_SUBSET_current_size
-! mapping from source number in current subset to real source number in all the subsets
- isource = isources_already_done + isource_in_this_subset
+ ! mapping from source number in current subset to real source number in all the subsets
+ isource = isources_already_done + isource_in_this_subset
-! loop on all the results to determine the best slice
- distmin = HUGEVAL
- do iprocloop = 0,NPROCTOT-1
- if(final_distance_source_all(isource_in_this_subset,iprocloop) < distmin) then
- distmin = final_distance_source_all(isource_in_this_subset,iprocloop)
- islice_selected_source(isource) = iprocloop
- ispec_selected_source(isource) = ispec_selected_source_all(isource_in_this_subset,iprocloop)
- xi_source(isource) = xi_source_all(isource_in_this_subset,iprocloop)
- eta_source(isource) = eta_source_all(isource_in_this_subset,iprocloop)
- gamma_source(isource) = gamma_source_all(isource_in_this_subset,iprocloop)
- x_found_source(isource_in_this_subset) = x_found_source_all(isource_in_this_subset,iprocloop)
- y_found_source(isource_in_this_subset) = y_found_source_all(isource_in_this_subset,iprocloop)
- z_found_source(isource_in_this_subset) = z_found_source_all(isource_in_this_subset,iprocloop)
- endif
- enddo
- final_distance_source(isource) = distmin
+ ! loop on all the results to determine the best slice
+ distmin = HUGEVAL
+ do iprocloop = 0,NPROCTOT-1
+ if(final_distance_source_all(isource_in_this_subset,iprocloop) < distmin) then
+ distmin = final_distance_source_all(isource_in_this_subset,iprocloop)
+ islice_selected_source(isource) = iprocloop
+ ispec_selected_source(isource) = ispec_selected_source_all(isource_in_this_subset,iprocloop)
+ xi_source(isource) = xi_source_all(isource_in_this_subset,iprocloop)
+ eta_source(isource) = eta_source_all(isource_in_this_subset,iprocloop)
+ gamma_source(isource) = gamma_source_all(isource_in_this_subset,iprocloop)
+ x_found_source(isource_in_this_subset) = x_found_source_all(isource_in_this_subset,iprocloop)
+ y_found_source(isource_in_this_subset) = y_found_source_all(isource_in_this_subset,iprocloop)
+ z_found_source(isource_in_this_subset) = z_found_source_all(isource_in_this_subset,iprocloop)
+ endif
+ enddo
+ final_distance_source(isource) = distmin
- write(IMAIN,*)
- write(IMAIN,*) '*************************************'
- write(IMAIN,*) ' locating source ',isource
- write(IMAIN,*) '*************************************'
- write(IMAIN,*)
- write(IMAIN,*) 'source located in slice ',islice_selected_source(isource_in_this_subset)
- write(IMAIN,*) ' in element ',ispec_selected_source(isource_in_this_subset)
- write(IMAIN,*)
- ! different output for force point sources
- if(USE_FORCE_POINT_SOURCE) then
- write(IMAIN,*) ' i index of source in that element: ',nint(xi_source(isource))
- write(IMAIN,*) ' j index of source in that element: ',nint(eta_source(isource))
- write(IMAIN,*) ' k index of source in that element: ',nint(gamma_source(isource))
write(IMAIN,*)
- write(IMAIN,*) ' component direction: ',COMPONENT_FORCE_SOURCE
+ write(IMAIN,*) '*************************************'
+ write(IMAIN,*) ' locating source ',isource
+ write(IMAIN,*) '*************************************'
write(IMAIN,*)
- write(IMAIN,*) ' nu1 = ',nu_source(1,:,isource)
- write(IMAIN,*) ' nu2 = ',nu_source(2,:,isource)
- write(IMAIN,*) ' nu3 = ',nu_source(3,:,isource)
+ write(IMAIN,*) 'source located in slice ',islice_selected_source(isource_in_this_subset)
+ write(IMAIN,*) ' in element ',ispec_selected_source(isource_in_this_subset)
write(IMAIN,*)
- write(IMAIN,*) ' at (x,y,z) coordinates = ',x_found_source(isource_in_this_subset),&
- y_found_source(isource_in_this_subset),z_found_source(isource_in_this_subset)
+ ! different output for force point sources
+ if(USE_FORCE_POINT_SOURCE) then
+ write(IMAIN,*) ' i index of source in that element: ',nint(xi_source(isource))
+ write(IMAIN,*) ' j index of source in that element: ',nint(eta_source(isource))
+ write(IMAIN,*) ' k index of source in that element: ',nint(gamma_source(isource))
+ write(IMAIN,*)
+ write(IMAIN,*) ' component direction: ',COMPONENT_FORCE_SOURCE
+ write(IMAIN,*)
+ write(IMAIN,*) ' nu1 = ',nu_source(1,:,isource)
+ write(IMAIN,*) ' nu2 = ',nu_source(2,:,isource)
+ write(IMAIN,*) ' nu3 = ',nu_source(3,:,isource)
+ write(IMAIN,*)
+ write(IMAIN,*) ' at (x,y,z) coordinates = ',x_found_source(isource_in_this_subset),&
+ y_found_source(isource_in_this_subset),z_found_source(isource_in_this_subset)
- ! prints frequency content for point forces
- f0 = hdur(isource)
- t0_ricker = 1.2d0/f0
- write(IMAIN,*) ' using a source of dominant frequency ',f0
- write(IMAIN,*) ' lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
- write(IMAIN,*) ' lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
- write(IMAIN,*) ' t0_ricker = ',t0_ricker,'tshift_cmt = ',tshift_cmt(isource)
- write(IMAIN,*)
- write(IMAIN,*) ' half duration -> frequency: ',hdur(isource),' seconds**(-1)'
- else
- write(IMAIN,*) ' xi coordinate of source in that element: ',xi_source(isource)
- write(IMAIN,*) ' eta coordinate of source in that element: ',eta_source(isource)
- write(IMAIN,*) 'gamma coordinate of source in that element: ',gamma_source(isource)
- ! add message if source is a Heaviside
- if(hdur(isource) <= 5.*DT) then
+ ! prints frequency content for point forces
+ f0 = hdur(isource)
+ t0_ricker = 1.2d0/f0
+ write(IMAIN,*) ' using a source of dominant frequency ',f0
+ write(IMAIN,*) ' lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ write(IMAIN,*) ' lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ write(IMAIN,*) ' t0_ricker = ',t0_ricker,'tshift_cmt = ',tshift_cmt(isource)
write(IMAIN,*)
- write(IMAIN,*) 'Source time function is a Heaviside, convolve later'
+ write(IMAIN,*) ' half duration -> frequency: ',hdur(isource),' seconds**(-1)'
+ else
+ write(IMAIN,*) ' xi coordinate of source in that element: ',xi_source(isource)
+ write(IMAIN,*) ' eta coordinate of source in that element: ',eta_source(isource)
+ write(IMAIN,*) 'gamma coordinate of source in that element: ',gamma_source(isource)
+ ! add message if source is a Heaviside
+ if(hdur(isource) <= 5.*DT) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'Source time function is a Heaviside, convolve later'
+ write(IMAIN,*)
+ endif
write(IMAIN,*)
+ write(IMAIN,*) ' half duration: ',hdur(isource),' seconds'
endif
- write(IMAIN,*)
- write(IMAIN,*) ' half duration: ',hdur(isource),' seconds'
- endif
- write(IMAIN,*) ' time shift: ',tshift_cmt(isource),' seconds'
+ write(IMAIN,*) ' time shift: ',tshift_cmt(isource),' seconds'
-! get latitude, longitude and depth of the source that will be used
- call xyz_2_rthetaphi_dble(x_found_source(isource_in_this_subset),y_found_source(isource_in_this_subset), &
+ ! writes out actual source position to vtk file
+ write(IOVTK,*) sngl(x_found_source(isource_in_this_subset)), &
+ sngl(y_found_source(isource_in_this_subset)), &
+ sngl(z_found_source(isource_in_this_subset))
+
+ ! get latitude, longitude and depth of the source that will be used
+ call xyz_2_rthetaphi_dble(x_found_source(isource_in_this_subset),y_found_source(isource_in_this_subset), &
z_found_source(isource_in_this_subset),r_found_source,theta_source(isource),phi_source(isource))
- call reduce(theta_source(isource),phi_source(isource))
+ call reduce(theta_source(isource),phi_source(isource))
-! convert geocentric to geographic colatitude
- colat_source = PI/2.0d0 &
+ ! convert geocentric to geographic colatitude
+ colat_source = PI/2.0d0 &
- datan(1.006760466d0*dcos(theta_source(isource))/dmax1(TINYVAL,dsin(theta_source(isource))))
- if(phi_source(isource)>PI) phi_source(isource)=phi_source(isource)-TWO_PI
+ if(phi_source(isource)>PI) phi_source(isource)=phi_source(isource)-TWO_PI
- write(IMAIN,*)
- write(IMAIN,*) 'original (requested) position of the source:'
- write(IMAIN,*)
- write(IMAIN,*) ' latitude: ',lat(isource)
- write(IMAIN,*) ' longitude: ',long(isource)
- write(IMAIN,*) ' depth: ',depth(isource),' km'
- write(IMAIN,*)
+ write(IMAIN,*)
+ write(IMAIN,*) 'original (requested) position of the source:'
+ write(IMAIN,*)
+ write(IMAIN,*) ' latitude: ',lat(isource)
+ write(IMAIN,*) ' longitude: ',long(isource)
+ write(IMAIN,*) ' depth: ',depth(isource),' km'
+ write(IMAIN,*)
-! compute real position of the source
- write(IMAIN,*) 'position of the source that will be used:'
- write(IMAIN,*)
- write(IMAIN,*) ' latitude: ',(PI/2.0d0-colat_source)*180.0d0/PI
- write(IMAIN,*) ' longitude: ',phi_source(isource)*180.0d0/PI
- write(IMAIN,*) ' depth: ',(r0-r_found_source)*R_EARTH/1000.0d0,' km'
- write(IMAIN,*)
+ ! compute real position of the source
+ write(IMAIN,*) 'position of the source that will be used:'
+ write(IMAIN,*)
+ write(IMAIN,*) ' latitude: ',(PI/2.0d0-colat_source)*180.0d0/PI
+ write(IMAIN,*) ' longitude: ',phi_source(isource)*180.0d0/PI
+ write(IMAIN,*) ' depth: ',(r0-r_found_source)*R_EARTH/1000.0d0,' km'
+ write(IMAIN,*)
-! display error in location estimate
- write(IMAIN,*) 'error in location of the source: ',sngl(final_distance_source(isource)),' km'
+ ! display error in location estimate
+ write(IMAIN,*) 'error in location of the source: ',sngl(final_distance_source(isource)),' km'
-! add warning if estimate is poor
-! (usually means source outside the mesh given by the user)
- if(final_distance_source(isource) > 50.d0) then
- write(IMAIN,*)
- write(IMAIN,*) '*****************************************************'
- write(IMAIN,*) '*****************************************************'
- write(IMAIN,*) '***** WARNING: source location estimate is poor *****'
- write(IMAIN,*) '*****************************************************'
- write(IMAIN,*) '*****************************************************'
- endif
+ ! add warning if estimate is poor
+ ! (usually means source outside the mesh given by the user)
+ if(final_distance_source(isource) > 50.d0) then
+ write(IMAIN,*)
+ write(IMAIN,*) '*****************************************************'
+ write(IMAIN,*) '*****************************************************'
+ write(IMAIN,*) '***** WARNING: source location estimate is poor *****'
+ write(IMAIN,*) '*****************************************************'
+ write(IMAIN,*) '*****************************************************'
+ endif
-! print source time function and spectrum
- if(PRINT_SOURCE_TIME_FUNCTION) then
+ ! print source time function and spectrum
+ if(PRINT_SOURCE_TIME_FUNCTION) then
- write(IMAIN,*)
- write(IMAIN,*) 'printing the source-time function'
+ write(IMAIN,*)
+ write(IMAIN,*) 'printing the source-time function'
- ! print the source-time function
- if(NSOURCES == 1) then
- plot_file = '/plot_source_time_function.txt'
- else
- if(isource < 10) then
- write(plot_file,"('/plot_source_time_function',i1,'.txt')") isource
- elseif(isource < 100) then
- write(plot_file,"('/plot_source_time_function',i2,'.txt')") isource
+ ! print the source-time function
+ if(NSOURCES == 1) then
+ plot_file = '/plot_source_time_function.txt'
else
- write(plot_file,"('/plot_source_time_function',i3,'.txt')") isource
+ if(isource < 10) then
+ write(plot_file,"('/plot_source_time_function',i1,'.txt')") isource
+ elseif(isource < 100) then
+ write(plot_file,"('/plot_source_time_function',i2,'.txt')") isource
+ else
+ write(plot_file,"('/plot_source_time_function',i3,'.txt')") isource
+ endif
endif
- endif
- open(unit=27,file=trim(OUTPUT_FILES)//plot_file,status='unknown')
+ open(unit=27,file=trim(OUTPUT_FILES)//plot_file,status='unknown')
- scalar_moment = 0.
- do i = 1,6
- scalar_moment = scalar_moment + moment_tensor(i,isource)**2
- enddo
- scalar_moment = dsqrt(scalar_moment/2.)
+ scalar_moment = 0.
+ do i = 1,6
+ scalar_moment = scalar_moment + moment_tensor(i,isource)**2
+ enddo
+ scalar_moment = dsqrt(scalar_moment/2.)
- ! define t0 as the earliest start time
- ! note: this calculation here is only used for outputting the plot_source_time_function file
- ! (see setup_sources_receivers.f90)
- t0 = - 1.5d0*minval( tshift_cmt(:) - hdur(:) )
- if( USE_FORCE_POINT_SOURCE ) t0 = - 1.2d0 * minval(tshift_cmt(:) - 1.0d0/hdur(:))
- t_cmt_used(:) = t_cmt_used(:)
- if( USER_T0 > 0.d0 ) then
- if( t0 <= USER_T0 + min_tshift_cmt_original ) then
- t_cmt_used(:) = tshift_cmt(:) + min_tshift_cmt_original
- t0 = USER_T0
+ ! define t0 as the earliest start time
+ ! note: this calculation here is only used for outputting the plot_source_time_function file
+ ! (see setup_sources_receivers.f90)
+ t0 = - 1.5d0*minval( tshift_cmt(:) - hdur(:) )
+ if( USE_FORCE_POINT_SOURCE ) t0 = - 1.2d0 * minval(tshift_cmt(:) - 1.0d0/hdur(:))
+ t_cmt_used(:) = t_cmt_used(:)
+ if( USER_T0 > 0.d0 ) then
+ if( t0 <= USER_T0 + min_tshift_cmt_original ) then
+ t_cmt_used(:) = tshift_cmt(:) + min_tshift_cmt_original
+ t0 = USER_T0
+ endif
endif
- endif
- ! convert the half duration for triangle STF to the one for gaussian STF
- ! note: this calculation here is only used for outputting the plot_source_time_function file
- ! (see setup_sources_receivers.f90)
- hdur_gaussian(:) = hdur(:)/SOURCE_DECAY_MIMIC_TRIANGLE
+ ! convert the half duration for triangle STF to the one for gaussian STF
+ ! note: this calculation here is only used for outputting the plot_source_time_function file
+ ! (see setup_sources_receivers.f90)
+ hdur_gaussian(:) = hdur(:)/SOURCE_DECAY_MIMIC_TRIANGLE
- ! writes out source time function to file
- do it=1,NSTEP
- time_source = dble(it-1)*DT-t0-t_cmt_used(isource)
- if( USE_FORCE_POINT_SOURCE ) then
- ! Ricker source time function
- f0 = hdur(isource)
- write(27,*) sngl(dble(it-1)*DT-t0), &
- sngl(FACTOR_FORCE_SOURCE*comp_source_time_function_rickr(time_source,f0))
- else
- ! Gaussian source time function
- write(27,*) sngl(dble(it-1)*DT-t0), &
- sngl(scalar_moment*comp_source_time_function(time_source,hdur_gaussian(isource)))
- endif
- enddo
- close(27)
+ ! writes out source time function to file
+ do it=1,NSTEP
+ time_source = dble(it-1)*DT-t0-t_cmt_used(isource)
+ if( USE_FORCE_POINT_SOURCE ) then
+ ! Ricker source time function
+ f0 = hdur(isource)
+ write(27,*) sngl(dble(it-1)*DT-t0), &
+ sngl(FACTOR_FORCE_SOURCE*comp_source_time_function_rickr(time_source,f0))
+ else
+ ! Gaussian source time function
+ write(27,*) sngl(dble(it-1)*DT-t0), &
+ sngl(scalar_moment*comp_source_time_function(time_source,hdur_gaussian(isource)))
+ endif
+ enddo
+ close(27)
- write(IMAIN,*)
- write(IMAIN,*) 'printing the source spectrum'
+ write(IMAIN,*)
+ write(IMAIN,*) 'printing the source spectrum'
- ! print the spectrum of the derivative of the source from 0 to 1/8 Hz
- if(NSOURCES == 1) then
- plot_file = '/plot_source_spectrum.txt'
- else
- if(isource < 10) then
- write(plot_file,"('/plot_source_spectrum',i1,'.txt')") isource
- elseif(isource < 100) then
- write(plot_file,"('/plot_source_spectrum',i2,'.txt')") isource
+ ! print the spectrum of the derivative of the source from 0 to 1/8 Hz
+ if(NSOURCES == 1) then
+ plot_file = '/plot_source_spectrum.txt'
else
- write(plot_file,"('/plot_source_spectrum',i3,'.txt')") isource
+ if(isource < 10) then
+ write(plot_file,"('/plot_source_spectrum',i1,'.txt')") isource
+ elseif(isource < 100) then
+ write(plot_file,"('/plot_source_spectrum',i2,'.txt')") isource
+ else
+ write(plot_file,"('/plot_source_spectrum',i3,'.txt')") isource
+ endif
endif
- endif
- open(unit=27,file=trim(OUTPUT_FILES)//plot_file,status='unknown')
+ open(unit=27,file=trim(OUTPUT_FILES)//plot_file,status='unknown')
- do iom=1,NSAMP_PLOT_SOURCE
- om=TWO_PI*(1.0d0/8.0d0)*(iom-1)/dble(NSAMP_PLOT_SOURCE-1)
- write(27,*) sngl(om/TWO_PI), &
- sngl(scalar_moment*om*comp_source_spectrum(om,hdur(isource)))
- enddo
- close(27)
+ do iom=1,NSAMP_PLOT_SOURCE
+ om=TWO_PI*(1.0d0/8.0d0)*(iom-1)/dble(NSAMP_PLOT_SOURCE-1)
+ write(27,*) sngl(om/TWO_PI), &
+ sngl(scalar_moment*om*comp_source_spectrum(om,hdur(isource)))
+ enddo
+ close(27)
- endif !PRINT_SOURCE_TIME_FUNCTION
+ endif !PRINT_SOURCE_TIME_FUNCTION
- enddo ! end of loop on all the sources within current source subset
+ enddo ! end of loop on all the sources within current source subset
endif ! end of section executed by main process only
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90 2011-04-02 04:05:46 UTC (rev 18164)
@@ -280,9 +280,12 @@
! save/read the surface movie using the same c routine as we do for absorbing boundaries (file ID is 9)
reclen=CUSTOM_REAL*NDIM*NGLLX*NGLLY*NSPEC_TOP*NSTEP
write(outputname,"('/proc',i6.6,'_surface_movie')") myrank
- if (NOISE_TOMOGRAPHY==1) call open_file_abs_w(9,trim(LOCAL_PATH)//trim(outputname),len_trim(trim(LOCAL_PATH)//trim(outputname)),reclen)
- if (NOISE_TOMOGRAPHY==2) call open_file_abs_r(9,trim(LOCAL_PATH)//trim(outputname),len_trim(trim(LOCAL_PATH)//trim(outputname)),reclen)
- if (NOISE_TOMOGRAPHY==3) call open_file_abs_r(9,trim(LOCAL_PATH)//trim(outputname),len_trim(trim(LOCAL_PATH)//trim(outputname)),reclen)
+ if (NOISE_TOMOGRAPHY==1) call open_file_abs_w(9,trim(LOCAL_PATH)//trim(outputname), &
+ len_trim(trim(LOCAL_PATH)//trim(outputname)),reclen)
+ if (NOISE_TOMOGRAPHY==2) call open_file_abs_r(9,trim(LOCAL_PATH)//trim(outputname), &
+ len_trim(trim(LOCAL_PATH)//trim(outputname)),reclen)
+ if (NOISE_TOMOGRAPHY==3) call open_file_abs_r(9,trim(LOCAL_PATH)//trim(outputname), &
+ len_trim(trim(LOCAL_PATH)//trim(outputname)),reclen)
endif
end subroutine check_parameters_noise
@@ -520,15 +523,19 @@
character(len=150) :: LOCAL_PATH
! output parameters
! local parameters
- integer :: ipoin,ispec2D,ispec,i,j,k,iglob
+ integer :: ispec2D,ispec,i,j,k,iglob
real(kind=CUSTOM_REAL), dimension(nmovie_points) :: &
store_val_x,store_val_y,store_val_z, &
store_val_ux,store_val_uy,store_val_uz
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: SURFACE_MOVIE
- character(len=150) :: outputname
+ !integer :: ipoin
+ !character(len=150) :: outputname
+ integer :: idummy
+ real(kind=CUSTOM_REAL) :: rdummy
+ character :: cdummy
+
+ allocate(SURFACE_MOVIE(NDIM,NGLLX,NGLLY,nspec_top))
-
- allocate(SURFACE_MOVIE(NDIM,NGLLX,NGLLY,nspec_top))
! get coordinates of surface mesh and surface displacement
do ispec2D = 1, nspec_top ! NSPEC2D_TOP(IREGION_CRUST_MANTLE)
ispec = ibelm_top_crust_mantle(ispec2D)
@@ -547,6 +554,20 @@
! change LOCAL_PATH specified in "DATA/Par_file"
call write_abs(9,SURFACE_MOVIE,CUSTOM_REAL*NDIM*NGLLX*NGLLY*nspec_top,it)
deallocate(SURFACE_MOVIE)
+
+ ! just to avoid compiler warnings
+ idummy = myrank
+ rdummy = xstore_crust_mantle(1)
+ rdummy = ystore_crust_mantle(1)
+ rdummy = zstore_crust_mantle(1)
+ rdummy = store_val_x(1)
+ rdummy = store_val_y(1)
+ rdummy = store_val_z(1)
+ rdummy = store_val_ux(1)
+ rdummy = store_val_uy(1)
+ rdummy = store_val_uz(1)
+ cdummy = LOCAL_PATH(1)
+
end subroutine noise_save_surface_movie
! =============================================================================================================
@@ -656,13 +677,15 @@
character(len=150) :: LOCAL_PATH
! output parameters
! local parameters
- integer :: ipoin,ispec2D,ispec,i,j,k,iglob,ios
+ integer :: ipoin,ispec2D,ispec,i,j,k,iglob !,ios
real(kind=CUSTOM_REAL), dimension(nmovie_points) :: store_val_ux,store_val_uy,store_val_uz
real(kind=CUSTOM_REAL) :: eta
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: SURFACE_MOVIE
- character(len=150) :: outputname
+ !character(len=150) :: outputname
+ integer :: idummy
+ real(kind=CUSTOM_REAL) :: rdummy
+ character :: cdummy
-
allocate(SURFACE_MOVIE(NDIM,NGLLX,NGLLY,nspec_top))
! read surface movie
call read_abs(9,SURFACE_MOVIE,CUSTOM_REAL*NDIM*NGLLX*NGLLY*nspec_top,it)
@@ -695,6 +718,14 @@
enddo
deallocate(SURFACE_MOVIE)
+
+ ! just to avoid compiler warnings
+ idummy = myrank
+ rdummy = store_val_ux(1)
+ rdummy = store_val_uy(1)
+ rdummy = store_val_uz(1)
+ cdummy = LOCAL_PATH(1)
+
end subroutine noise_read_add_surface_movie
! =============================================================================================================
@@ -798,13 +829,15 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
Sigma_kl_crust_mantle
! local parameters
- integer :: i,j,k,ispec,iglob,ipoin,ispec2D,ios
+ integer :: i,j,k,ispec,iglob,ipoin,ispec2D !,ios
real(kind=CUSTOM_REAL) :: eta
real(kind=CUSTOM_REAL), dimension(nmovie_points) :: store_val_ux,store_val_uy,store_val_uz
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: SURFACE_MOVIE
- character(len=150) :: outputname
-
-
+ !character(len=150) :: outputname
+ integer :: idummy
+ real(kind=CUSTOM_REAL) :: rdummy
+ character :: cdummy
+
allocate(SURFACE_MOVIE(NDIM,NGLLX,NGLLY,nspec_top))
! read surface movie, needed for Sigma_kl_crust_mantle
call read_abs(9,SURFACE_MOVIE,CUSTOM_REAL*NDIM*NGLLX*NGLLY*nspec_top,it)
@@ -836,6 +869,14 @@
enddo
+ ! just to avoid compiler warnings
+ idummy = myrank
+ rdummy = store_val_ux(1)
+ rdummy = store_val_uy(1)
+ rdummy = store_val_uz(1)
+ cdummy = LOCAL_PATH(1)
+
+
end subroutine compute_kernels_strength_noise
! =============================================================================================================
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90 2011-04-02 04:05:46 UTC (rev 18164)
@@ -333,12 +333,29 @@
close(IOVTK)
! we should know NSOURCES+nrec at this point...
+ ! creates source/receiver location file
filename = trim(OUTPUT_FILES)//'/sr_tmp.vtk'
filename_new = trim(OUTPUT_FILES)//'/sr.vtk'
- write(system_command,"('sed -e ',a1,'s/POINTS.*/POINTS',i6,' float/',a1,' < ',a,' > ',a)") &
+ write(system_command, &
+ "('sed -e ',a1,'s/POINTS.*/POINTS',i6,' float/',a1,' < ',a,' > ',a)")&
"'",NSOURCES + nrec,"'",trim(filename),trim(filename_new)
call system(system_command)
+ ! only extract receiver locations and remove temporary file
+ filename_new = trim(OUTPUT_FILES)//'/receiver.vtk'
+ write(system_command, &
+ "('awk ',a1,'{if(NR<5) print $0;if(NR==6)print ',a1,'POINTS',i6,' float',a1,';if(NR>5+',i6,')print $0}',a1,' < ',a,' > ',a)")&
+ "'",'"',nrec,'"',NSOURCES,"'",trim(filename),trim(filename_new)
+ call system(system_command)
+
+ ! only extract source locations and remove temporary file
+ filename_new = trim(OUTPUT_FILES)//'/source.vtk'
+ write(system_command, &
+ "('awk ',a1,'{if(NR< 6 + ',i6,') print $0}END{print}',a1,' < ',a,' > ',a,'; rm -f ',a)")&
+ "'",NSOURCES,"'",trim(filename),trim(filename_new),trim(filename)
+ call system(system_command)
+
+
write(IMAIN,*)
write(IMAIN,*) 'Total number of samples for seismograms = ',NSTEP
write(IMAIN,*)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90 2011-04-02 04:05:46 UTC (rev 18164)
@@ -1066,11 +1066,13 @@
npoin2D_max_all_CM_IC = max(maxval(npoin2D_xi_crust_mantle(:) + npoin2D_xi_inner_core(:)), &
maxval(npoin2D_eta_crust_mantle(:) + npoin2D_eta_inner_core(:)))
- allocate(buffer_send_faces(NDIM,npoin2D_max_all_CM_IC,NUMFACES_SHARED))
- allocate(buffer_received_faces(NDIM,npoin2D_max_all_CM_IC,NUMFACES_SHARED))
+ allocate(buffer_send_faces(NDIM,npoin2D_max_all_CM_IC,NUMFACES_SHARED), &
+ buffer_received_faces(NDIM,npoin2D_max_all_CM_IC,NUMFACES_SHARED),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating mpi buffer')
- allocate(b_buffer_send_faces(NDIM,npoin2D_max_all_CM_IC,NUMFACES_SHARED))
- allocate(b_buffer_received_faces(NDIM,npoin2D_max_all_CM_IC,NUMFACES_SHARED))
+ allocate(b_buffer_send_faces(NDIM,npoin2D_max_all_CM_IC,NUMFACES_SHARED), &
+ b_buffer_received_faces(NDIM,npoin2D_max_all_CM_IC,NUMFACES_SHARED),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating mpi b_buffer')
call fix_non_blocking_slices(is_on_a_slice_edge_crust_mantle,iboolright_xi_crust_mantle, &
iboolleft_xi_crust_mantle,iboolright_eta_crust_mantle,iboolleft_eta_crust_mantle, &
@@ -1096,7 +1098,8 @@
else
nabs_xmin_cm = 1
endif
- allocate(absorb_xmin_crust_mantle5(NDIM,NGLLY,NGLLZ,nabs_xmin_cm,8))
+ allocate(absorb_xmin_crust_mantle5(NDIM,NGLLY,NGLLZ,nabs_xmin_cm,8),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmin')
if (nspec2D_xmax_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
@@ -1104,7 +1107,8 @@
else
nabs_xmax_cm = 1
endif
- allocate(absorb_xmax_crust_mantle5(NDIM,NGLLY,NGLLZ,nabs_xmax_cm,8))
+ allocate(absorb_xmax_crust_mantle5(NDIM,NGLLY,NGLLZ,nabs_xmax_cm,8),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmax')
if (nspec2D_ymin_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
@@ -1112,7 +1116,8 @@
else
nabs_ymin_cm = 1
endif
- allocate(absorb_ymin_crust_mantle5(NDIM,NGLLX,NGLLZ,nabs_ymin_cm,8))
+ allocate(absorb_ymin_crust_mantle5(NDIM,NGLLX,NGLLZ,nabs_ymin_cm,8),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymin')
if (nspec2D_ymax_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
@@ -1120,7 +1125,8 @@
else
nabs_ymax_cm = 1
endif
- allocate(absorb_ymax_crust_mantle5(NDIM,NGLLX,NGLLZ,nabs_ymax_cm,8))
+ allocate(absorb_ymax_crust_mantle5(NDIM,NGLLX,NGLLZ,nabs_ymax_cm,8),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymax')
! outer_core
if (nspec2D_xmin_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
@@ -1129,7 +1135,8 @@
else
nabs_xmin_oc = 1
endif
- allocate(absorb_xmin_outer_core(NGLLY,NGLLZ,nabs_xmin_oc))
+ allocate(absorb_xmin_outer_core(NGLLY,NGLLZ,nabs_xmin_oc),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmin')
if (nspec2D_xmax_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
@@ -1137,7 +1144,8 @@
else
nabs_xmax_oc = 1
endif
- allocate(absorb_xmax_outer_core(NGLLY,NGLLZ,nabs_xmax_oc))
+ allocate(absorb_xmax_outer_core(NGLLY,NGLLZ,nabs_xmax_oc),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmax')
if (nspec2D_ymin_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
@@ -1145,7 +1153,8 @@
else
nabs_ymin_oc = 1
endif
- allocate(absorb_ymin_outer_core(NGLLX,NGLLZ,nabs_ymin_oc))
+ allocate(absorb_ymin_outer_core(NGLLX,NGLLZ,nabs_ymin_oc),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymin')
if (nspec2D_ymax_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
@@ -1153,7 +1162,8 @@
else
nabs_ymax_oc = 1
endif
- allocate(absorb_ymax_outer_core(NGLLX,NGLLZ,nabs_ymax_oc))
+ allocate(absorb_ymax_outer_core(NGLLX,NGLLZ,nabs_ymax_oc),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymax')
if (NSPEC2D_BOTTOM(IREGION_OUTER_CORE) > 0 .and. &
(SIMULATION_TYPE == 3 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
@@ -1161,7 +1171,8 @@
else
nabs_zmin_oc = 1
endif
- allocate(absorb_zmin_outer_core(NGLLX,NGLLY,nabs_zmin_oc))
+ allocate(absorb_zmin_outer_core(NGLLX,NGLLY,nabs_zmin_oc),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb zmin')
! read arrays for Stacey conditions
call read_mesh_databases_stacey(myrank, &
@@ -1193,37 +1204,39 @@
! source and receivers
! allocate arrays for source
- allocate(islice_selected_source(NSOURCES))
- allocate(ispec_selected_source(NSOURCES))
- allocate(Mxx(NSOURCES))
- allocate(Myy(NSOURCES))
- allocate(Mzz(NSOURCES))
- allocate(Mxy(NSOURCES))
- allocate(Mxz(NSOURCES))
- allocate(Myz(NSOURCES))
- allocate(xi_source(NSOURCES))
- allocate(eta_source(NSOURCES))
- allocate(gamma_source(NSOURCES))
- allocate(tshift_cmt(NSOURCES))
- allocate(hdur(NSOURCES))
- allocate(hdur_gaussian(NSOURCES))
- allocate(theta_source(NSOURCES))
- allocate(phi_source(NSOURCES))
- allocate(nu_source(NDIM,NDIM,NSOURCES))
+ allocate(islice_selected_source(NSOURCES), &
+ ispec_selected_source(NSOURCES), &
+ Mxx(NSOURCES), &
+ Myy(NSOURCES), &
+ Mzz(NSOURCES), &
+ Mxy(NSOURCES), &
+ Mxz(NSOURCES), &
+ Myz(NSOURCES), &
+ xi_source(NSOURCES), &
+ eta_source(NSOURCES), &
+ gamma_source(NSOURCES), &
+ tshift_cmt(NSOURCES), &
+ hdur(NSOURCES), &
+ hdur_gaussian(NSOURCES), &
+ theta_source(NSOURCES), &
+ phi_source(NSOURCES), &
+ nu_source(NDIM,NDIM,NSOURCES),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating source arrays')
! allocate memory for receiver arrays
- allocate(islice_selected_rec(nrec))
- allocate(ispec_selected_rec(nrec))
- allocate(xi_receiver(nrec))
- allocate(eta_receiver(nrec))
- allocate(gamma_receiver(nrec))
- allocate(station_name(nrec))
- allocate(network_name(nrec))
- allocate(stlat(nrec))
- allocate(stlon(nrec))
- allocate(stele(nrec))
- allocate(stbur(nrec))
- allocate(nu(NDIM,NDIM,nrec))
+ allocate(islice_selected_rec(nrec), &
+ ispec_selected_rec(nrec), &
+ xi_receiver(nrec), &
+ eta_receiver(nrec), &
+ gamma_receiver(nrec), &
+ station_name(nrec), &
+ network_name(nrec), &
+ stlat(nrec), &
+ stlon(nrec), &
+ stele(nrec), &
+ stbur(nrec), &
+ nu(NDIM,NDIM,nrec),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating receiver arrays')
! locates sources and receivers
call setup_sources_receivers(NSOURCES,myrank,ibool_crust_mantle, &
@@ -1243,7 +1256,8 @@
! allocates source arrays
if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then
- allocate(sourcearrays(NDIM,NGLLX,NGLLY,NGLLZ,NSOURCES))
+ allocate(sourcearrays(NDIM,NGLLX,NGLLY,NGLLZ,NSOURCES),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating sourcearrays')
! stores source arrays
call setup_sources_receivers_srcarr(NSOURCES,myrank, &
@@ -1259,7 +1273,9 @@
if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
NSTEP_SUB_ADJ = ceiling( dble(NSTEP)/dble(NTSTEP_BETWEEN_READ_ADJSRC) )
- allocate(iadj_vec(NSTEP))
+ allocate(iadj_vec(NSTEP),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating iadj_vec')
+
! initializes iadj_vec
do it=1,NSTEP
iadj_vec(it) = NSTEP-it+1 ! default is for reversing entire record
@@ -1267,12 +1283,14 @@
if(nadj_rec_local > 0) then
! allocate adjoint source arrays
- allocate(adj_sourcearrays(NDIM,NGLLX,NGLLY,NGLLZ,nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC))
+ allocate(adj_sourcearrays(NDIM,NGLLX,NGLLY,NGLLZ,nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating adjoint sourcearrays')
adj_sourcearrays = 0._CUSTOM_REAL
! allocate indexing arrays
- allocate(iadjsrc(NSTEP_SUB_ADJ,2))
- allocate(iadjsrc_len(NSTEP_SUB_ADJ))
+ allocate(iadjsrc(NSTEP_SUB_ADJ,2), &
+ iadjsrc_len(NSTEP_SUB_ADJ),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating adjoint indexing arrays')
! initializes iadjsrc, iadjsrc_len and iadj_vec
call setup_sources_receivers_adjindx(NSTEP,NSTEP_SUB_ADJ, &
NTSTEP_BETWEEN_READ_ADJSRC, &
@@ -1283,20 +1301,23 @@
! allocates receiver interpolators
if (nrec_local > 0) then
! allocate Lagrange interpolators for receivers
- allocate(hxir_store(nrec_local,NGLLX))
- allocate(hetar_store(nrec_local,NGLLY))
- allocate(hgammar_store(nrec_local,NGLLZ))
+ allocate(hxir_store(nrec_local,NGLLX), &
+ hetar_store(nrec_local,NGLLY), &
+ hgammar_store(nrec_local,NGLLZ),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating receiver interpolators')
! define local to global receiver numbering mapping
- allocate(number_receiver_global(nrec_local))
+ allocate(number_receiver_global(nrec_local),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating global receiver numbering')
! define and store Lagrange interpolators at all the receivers
if (SIMULATION_TYPE == 2) then
nadj_hprec_local = nrec_local
else
nadj_hprec_local = 1
endif
- allocate(hpxir_store(nadj_hprec_local,NGLLX))
- allocate(hpetar_store(nadj_hprec_local,NGLLY))
- allocate(hpgammar_store(nadj_hprec_local,NGLLZ))
+ allocate(hpxir_store(nadj_hprec_local,NGLLX), &
+ hpetar_store(nadj_hprec_local,NGLLY), &
+ hpgammar_store(nadj_hprec_local,NGLLZ),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating derivative interpolators')
! stores interpolators for receiver positions
call setup_sources_receivers_intp(NSOURCES,myrank, &
@@ -1317,7 +1338,10 @@
allocate(seismograms(NDIM*NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
if(ier /= 0) stop 'error while allocating seismograms'
! allocate Frechet derivatives array
- allocate(moment_der(NDIM,NDIM,nrec_local),sloc_der(NDIM,nrec_local),stshift_der(nrec_local),shdur_der(nrec_local))
+ allocate(moment_der(NDIM,NDIM,nrec_local),sloc_der(NDIM,nrec_local), &
+ stshift_der(nrec_local),shdur_der(nrec_local),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating frechet derivatives arrays')
+
moment_der = 0._CUSTOM_REAL
sloc_der = 0._CUSTOM_REAL
stshift_der = 0._CUSTOM_REAL
@@ -1456,13 +1480,14 @@
endif
! allocate buffers for cube and slices
- allocate(sender_from_slices_to_cube(non_zero_nb_msgs_theor_in_cube))
- allocate(buffer_all_cube_from_slices(non_zero_nb_msgs_theor_in_cube,npoin2D_cube_from_slices,NDIM))
- allocate(b_buffer_all_cube_from_slices(non_zero_nb_msgs_theor_in_cube,npoin2D_cube_from_slices,NDIM))
- allocate(buffer_slices(npoin2D_cube_from_slices,NDIM))
- allocate(b_buffer_slices(npoin2D_cube_from_slices,NDIM))
- allocate(buffer_slices2(npoin2D_cube_from_slices,NDIM))
- allocate(ibool_central_cube(non_zero_nb_msgs_theor_in_cube,npoin2D_cube_from_slices))
+ allocate(sender_from_slices_to_cube(non_zero_nb_msgs_theor_in_cube), &
+ buffer_all_cube_from_slices(non_zero_nb_msgs_theor_in_cube,npoin2D_cube_from_slices,NDIM), &
+ b_buffer_all_cube_from_slices(non_zero_nb_msgs_theor_in_cube,npoin2D_cube_from_slices,NDIM), &
+ buffer_slices(npoin2D_cube_from_slices,NDIM), &
+ b_buffer_slices(npoin2D_cube_from_slices,NDIM), &
+ buffer_slices2(npoin2D_cube_from_slices,NDIM), &
+ ibool_central_cube(non_zero_nb_msgs_theor_in_cube,npoin2D_cube_from_slices),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating cube buffers')
! handles the communications with the central cube if it was included in the mesh
call prepare_timerun_centralcube(myrank,rmass_inner_core, &
@@ -1488,13 +1513,14 @@
! allocate fictitious buffers for cube and slices with a dummy size
! just to be able to use them as arguments in subroutine calls
- allocate(sender_from_slices_to_cube(1))
- allocate(buffer_all_cube_from_slices(1,1,1))
- allocate(b_buffer_all_cube_from_slices(1,1,1))
- allocate(buffer_slices(1,1))
- allocate(b_buffer_slices(1,1))
- allocate(buffer_slices2(1,1))
- allocate(ibool_central_cube(1,1))
+ allocate(sender_from_slices_to_cube(1), &
+ buffer_all_cube_from_slices(1,1,1), &
+ b_buffer_all_cube_from_slices(1,1,1), &
+ buffer_slices(1,1), &
+ b_buffer_slices(1,1), &
+ buffer_slices2(1,1), &
+ ibool_central_cube(1,1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating dummy buffers')
endif
@@ -1557,20 +1583,34 @@
nmovie_points = NGLLX * NGLLY * NSPEC2D_TOP(IREGION_CRUST_MANTLE)
NIT = 1
endif
- allocate(store_val_x(nmovie_points))
- allocate(store_val_y(nmovie_points))
- allocate(store_val_z(nmovie_points))
- allocate(store_val_ux(nmovie_points))
- allocate(store_val_uy(nmovie_points))
- allocate(store_val_uz(nmovie_points))
+ allocate(store_val_x(nmovie_points), &
+ store_val_y(nmovie_points), &
+ store_val_z(nmovie_points), &
+ store_val_ux(nmovie_points), &
+ store_val_uy(nmovie_points), &
+ store_val_uz(nmovie_points),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating movie surface arrays')
+
if (MOVIE_SURFACE) then ! those arrays are not neccessary for noise tomography, so only allocate them in MOVIE_SURFACE case
- allocate(store_val_x_all(nmovie_points,0:NPROCTOT_VAL-1))
- allocate(store_val_y_all(nmovie_points,0:NPROCTOT_VAL-1))
- allocate(store_val_z_all(nmovie_points,0:NPROCTOT_VAL-1))
- allocate(store_val_ux_all(nmovie_points,0:NPROCTOT_VAL-1))
- allocate(store_val_uy_all(nmovie_points,0:NPROCTOT_VAL-1))
- allocate(store_val_uz_all(nmovie_points,0:NPROCTOT_VAL-1))
+ allocate(store_val_x_all(nmovie_points,0:NPROCTOT_VAL-1), &
+ store_val_y_all(nmovie_points,0:NPROCTOT_VAL-1), &
+ store_val_z_all(nmovie_points,0:NPROCTOT_VAL-1), &
+ store_val_ux_all(nmovie_points,0:NPROCTOT_VAL-1), &
+ store_val_uy_all(nmovie_points,0:NPROCTOT_VAL-1), &
+ store_val_uz_all(nmovie_points,0:NPROCTOT_VAL-1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating movie surface all arrays')
endif
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'Movie surface:'
+ write(IMAIN,*) ' Writing to moviedata*** files in output directory'
+ if(MOVIE_VOLUME_TYPE == 5) then
+ write(IMAIN,*) ' movie output: displacement'
+ else
+ write(IMAIN,*) ' movie output: velocity'
+ endif
+ write(IMAIN,*) ' time steps every: ',NTSTEP_BETWEEN_FRAMES
+ endif
endif
@@ -1589,18 +1629,34 @@
MOVIE_COARSE,npoints_3dmovie,nspecel_3dmovie,num_ibool_3dmovie,mask_ibool,mask_3dmovie)
- allocate(nu_3dmovie(3,3,npoints_3dmovie))
+ allocate(nu_3dmovie(3,3,npoints_3dmovie),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating nu for 3d movie')
call write_movie_volume_mesh(npoints_3dmovie,prname,ibool_crust_mantle,xstore_crust_mantle, &
ystore_crust_mantle,zstore_crust_mantle, muvstore_crust_mantle_3dmovie, &
mask_3dmovie,mask_ibool,num_ibool_3dmovie,nu_3dmovie,MOVIE_COARSE)
if(myrank == 0) then
- write(IMAIN,*) 'Writing to movie3D files on local disk'
- write(IMAIN,*) 'depth(T,B):',MOVIE_TOP,MOVIE_BOTTOM
- write(IMAIN,*) 'lon(W,E) :',MOVIE_WEST,MOVIE_EAST
- write(IMAIN,*) 'lat(S,N) :',MOVIE_SOUTH,MOVIE_NORTH
- write(IMAIN,*) 'Starting at time step:',MOVIE_START, 'ending at:',MOVIE_STOP,'every: ',NTSTEP_BETWEEN_FRAMES
+ write(IMAIN,*)
+ write(IMAIN,*) 'Movie volume:'
+ write(IMAIN,*) ' Writing to movie3D*** files on local disk databases directory'
+ if(MOVIE_VOLUME_TYPE == 1) then
+ write(IMAIN,*) ' movie output: strain'
+ else if(MOVIE_VOLUME_TYPE == 2) then
+ write(IMAIN,*) ' movie output: time integral of strain'
+ else if(MOVIE_VOLUME_TYPE == 3) then
+ write(IMAIN,*) ' movie output: potency or integral of strain'
+ else if(MOVIE_VOLUME_TYPE == 4) then
+ write(IMAIN,*) ' movie output: divergence and curl'
+ else if(MOVIE_VOLUME_TYPE == 5) then
+ write(IMAIN,*) ' movie output: displacement'
+ else if(MOVIE_VOLUME_TYPE == 6) then
+ write(IMAIN,*) ' movie output: velocity'
+ endif
+ write(IMAIN,*) ' depth(T,B):',MOVIE_TOP,MOVIE_BOTTOM
+ write(IMAIN,*) ' lon(W,E) :',MOVIE_WEST,MOVIE_EAST
+ write(IMAIN,*) ' lat(S,N) :',MOVIE_SOUTH,MOVIE_NORTH
+ write(IMAIN,*) ' Starting at time step:',MOVIE_START, 'ending at:',MOVIE_STOP,'every: ',NTSTEP_BETWEEN_FRAMES
endif
endif ! MOVIE_VOLUME
@@ -1695,7 +1751,8 @@
! approximate hessian
if( APPROXIMATE_HESS_KL ) then
- allocate( hess_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT))
+ allocate( hess_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating hessian')
hess_kl_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
endif
@@ -1718,7 +1775,8 @@
else
nspec_beta_kl_outer_core = 1
endif
- allocate(beta_kl_outer_core(NGLLX,NGLLY,NGLLZ,nspec_beta_kl_outer_core))
+ allocate(beta_kl_outer_core(NGLLX,NGLLY,NGLLZ,nspec_beta_kl_outer_core),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating beta outercore')
beta_kl_outer_core(:,:,:,:) = 0._CUSTOM_REAL
endif
@@ -1774,11 +1832,13 @@
!<YANGL
! NOISE TOMOGRAPHY
if ( NOISE_TOMOGRAPHY /= 0 ) then
- allocate(noise_sourcearray(NDIM,NGLLX,NGLLY,NGLLZ,NSTEP))
- allocate(normal_x_noise(nmovie_points))
- allocate(normal_y_noise(nmovie_points))
- allocate(normal_z_noise(nmovie_points))
- allocate(mask_noise(nmovie_points))
+ allocate(noise_sourcearray(NDIM,NGLLX,NGLLY,NGLLZ,NSTEP), &
+ normal_x_noise(nmovie_points), &
+ normal_y_noise(nmovie_points), &
+ normal_z_noise(nmovie_points), &
+ mask_noise(nmovie_points),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating noise arrays')
+
noise_sourcearray(:,:,:,:,:) = 0._CUSTOM_REAL
normal_x_noise(:) = 0._CUSTOM_REAL
normal_y_noise(:) = 0._CUSTOM_REAL
@@ -4177,6 +4237,7 @@
if( mod(it,NTSTEP_BETWEEN_FRAMES) == 0) then
! save velocity here to avoid static offset on displacement for movies
call write_movie_surface(myrank,nmovie_points,scale_veloc,veloc_crust_mantle, &
+ scale_displ,displ_crust_mantle, &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
store_val_x,store_val_y,store_val_z, &
store_val_x_all,store_val_y_all,store_val_z_all, &
@@ -4184,7 +4245,7 @@
store_val_ux_all,store_val_uy_all,store_val_uz_all, &
ibelm_top_crust_mantle,ibool_crust_mantle, &
NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
- NIT,it,OUTPUT_FILES)
+ NIT,it,OUTPUT_FILES,MOVIE_VOLUME_TYPE)
endif
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_surface.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_surface.f90 2011-04-01 23:10:54 UTC (rev 18163)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_surface.f90 2011-04-02 04:05:46 UTC (rev 18164)
@@ -26,13 +26,14 @@
!=====================================================================
subroutine write_movie_surface(myrank,nmovie_points,scale_veloc,veloc_crust_mantle, &
+ scale_displ,displ_crust_mantle, &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
store_val_x,store_val_y,store_val_z, &
store_val_x_all,store_val_y_all,store_val_z_all, &
store_val_ux,store_val_uy,store_val_uz, &
store_val_ux_all,store_val_uy_all,store_val_uz_all, &
ibelm_top_crust_mantle,ibool_crust_mantle,nspec_top, &
- NIT,it,OUTPUT_FILES)
+ NIT,it,OUTPUT_FILES,MOVIE_VOLUME_TYPE)
implicit none
@@ -42,10 +43,10 @@
include "OUTPUT_FILES/values_from_mesher.h"
integer myrank,nmovie_points
- double precision :: scale_veloc
+ double precision :: scale_veloc,scale_displ
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
- veloc_crust_mantle
+ veloc_crust_mantle,displ_crust_mantle
real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle
@@ -64,6 +65,8 @@
integer nspec_top,NIT,it
character(len=150) OUTPUT_FILES
+ integer MOVIE_VOLUME_TYPE
+
! local parameters
character(len=150) :: outputname
integer :: ipoin,ispec2D,ispec,i,j,k,ier,iglob
@@ -90,9 +93,18 @@
store_val_x(ipoin) = xstore_crust_mantle(iglob)
store_val_y(ipoin) = ystore_crust_mantle(iglob)
store_val_z(ipoin) = zstore_crust_mantle(iglob)
- store_val_ux(ipoin) = veloc_crust_mantle(1,iglob)*scale_veloc
- store_val_uy(ipoin) = veloc_crust_mantle(2,iglob)*scale_veloc
- store_val_uz(ipoin) = veloc_crust_mantle(3,iglob)*scale_veloc
+ if(MOVIE_VOLUME_TYPE == 5) then
+ ! stores displacement
+ store_val_ux(ipoin) = displ_crust_mantle(1,iglob)*scale_displ
+ store_val_uy(ipoin) = displ_crust_mantle(2,iglob)*scale_displ
+ store_val_uz(ipoin) = displ_crust_mantle(3,iglob)*scale_displ
+ else
+ ! stores velocity
+ store_val_ux(ipoin) = veloc_crust_mantle(1,iglob)*scale_veloc
+ store_val_uy(ipoin) = veloc_crust_mantle(2,iglob)*scale_veloc
+ store_val_uz(ipoin) = veloc_crust_mantle(3,iglob)*scale_veloc
+ endif
+
enddo
enddo
More information about the CIG-COMMITS
mailing list