[cig-commits] [commit] devel: added write_output_SU.f90 (ab5df80)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Dec 10 12:05:27 PST 2014
Repository : https://github.com/geodynamics/specfem2d
On branch : devel
Link : https://github.com/geodynamics/specfem2d/compare/ee73d7f81d9311271516b8a722ca6f2f50ae8743...0509fc8b3f04a8436eaccda15f30b6cee82c654f
>---------------------------------------------------------------
commit ab5df80c477cc105d6e2114111884e4edeafaf59
Author: rmodrak <rmodrak at princeton.edu>
Date: Wed Dec 10 10:12:32 2014 -0500
added write_output_SU.f90
>---------------------------------------------------------------
ab5df80c477cc105d6e2114111884e4edeafaf59
src/specfem2D/Makefile.in | 5 +++
src/specfem2D/write_output_SU.f90 | 79 +++++++++++++++++++++++++++++++++++++
src/specfem2D/write_seismograms.F90 | 54 +------------------------
3 files changed, 85 insertions(+), 53 deletions(-)
diff --git a/src/specfem2D/Makefile.in b/src/specfem2D/Makefile.in
index 29163f0..5a65240 100644
--- a/src/specfem2D/Makefile.in
+++ b/src/specfem2D/Makefile.in
@@ -163,6 +163,7 @@ OBJS_SPECFEM2D = \
$O/setup_sources_receivers.o \
$O/sort_array_coordinates.o \
$O/update_displacement_scheme.o \
+ $O/write_output_SU.o \
$O/write_seismograms.o \
$O/specfem2D.o \
$O/write_jpeg_image.o \
@@ -291,6 +292,7 @@ $O/compute_curl_one_element.o: $O/specfem2D_par.o
$O/compute_pressure.o: $O/specfem2D_par.o
$O/compute_vector_field.o: $O/specfem2D_par.o
$O/create_color_image.o: $O/specfem2D_par.o
+$O/write_output_SU.o: $O/specfem2D_par.o
$O/write_seismograms.o: $O/specfem2D_par.o
$O/plotpost.o: $O/specfem2D_par.o
$O/prepare_color_image.o: $O/specfem2D_par.o
@@ -483,6 +485,9 @@ $O/sort_array_coordinates.o: ${S}/sort_array_coordinates.F90 ${SETUP}/constants.
$O/update_displacement_scheme.o: ${S}/update_displacement_scheme.F90 ${SETUP}/constants.h
${F90} $(DEF_FFLAGS) -c -o $O/update_displacement_scheme.o ${S}/update_displacement_scheme.F90
+$O/write_output_SU.o: ${S}/write_output_SU.f90 ${SETUP}/constants.h
+ ${F90} $(DEF_FFLAGS) -c -o $O/write_output_SU.o ${S}/write_output_SU.f90
+
$O/write_seismograms.o: ${S}/write_seismograms.F90 ${SETUP}/constants.h
${F90} $(DEF_FFLAGS) -c -o $O/write_seismograms.o ${S}/write_seismograms.F90
diff --git a/src/specfem2D/write_output_SU.f90 b/src/specfem2D/write_output_SU.f90
new file mode 100644
index 0000000..d815567
--- /dev/null
+++ b/src/specfem2D/write_output_SU.f90
@@ -0,0 +1,79 @@
+ subroutine write_output_SU(x_source,z_source,irec,buffer_binary,number_of_components)
+
+ use specfem_par, only : sisux,sisuz,siscurl,station_name,network_name, &
+ NSTEP,nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
+ NSTEP_BETWEEN_OUTPUT_SEISMOS,seismo_offset,seismo_current,p_sv, &
+ st_zval,SU_FORMAT,save_ASCII_seismograms, &
+ save_binary_seismograms_double,subsamp_seismos
+ integer :: deltat_int2
+ integer :: irec,isample,number_of_components
+
+! to write seismograms in single precision SEP and double precision binary
+! format
+ double precision, dimension(NSTEP_BETWEEN_OUTPUT_SEISMOS/subsamp_seismos,number_of_components) :: buffer_binary
+
+! scaling factor for Seismic Unix xsu dislay
+ double precision, parameter :: FACTORXSU = 1.d0
+
+
+ double precision :: x_source,z_source
+ integer(kind=2) :: header2(2)
+
+ print*, 'x_source', x_source
+ print*, 'z_source', z_source
+
+
+ if (seismo_offset==0) then
+
+ if (deltat*1.0d6 > 2**15) then
+ deltat_int2 = 0
+ else
+ deltat_int2 = NINT(deltat*1.0d6, kind=2) ! deltat (unit: 10^{-6} second)
+ endif
+
+ print*, 'made it here-1'
+
+ ! 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)=int(NSTEP, kind=2)
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+29) header2
+ header2(1)=deltat_int2
+ header2(2)=0 ! dummy
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+30) header2
+ 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)=int(NSTEP, kind=2)
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+29) header2
+ header2(1)=deltat_int2
+ header2(2)=0 ! dummy
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+30) header2
+ endif
+ endif
+
+ endif
+
+ ! the "60" in the following corresponds to 240 bytes header (note the reclength is 4 bytes)
+ do isample = 1, seismo_current
+ write(12,rec=irec*60+(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,1))
+ if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then
+ write(14,rec=irec*60+(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2))
+ endif
+ enddo
+
+ end subroutine write_output_SU
diff --git a/src/specfem2D/write_seismograms.F90 b/src/specfem2D/write_seismograms.F90
index 036686e..a5c5326 100644
--- a/src/specfem2D/write_seismograms.F90
+++ b/src/specfem2D/write_seismograms.F90
@@ -65,7 +65,6 @@
include "constants.h"
- integer :: deltat_int2
logical :: save_binary_seismograms
integer irec,length_station_name,length_network_name,iorientation,isample,number_of_components
@@ -82,10 +81,7 @@
integer :: irecloc
-!<SU_FORMAT
double precision :: x_source,z_source
- integer(kind=2) :: header2(2)
-!>SU_FORMAT
#ifdef USE_MPI
integer :: ierror
@@ -324,55 +320,7 @@
else ! if SU_FORMAT
- if (seismo_offset==0) then
-
- if (deltat*1.0d6 > 2**15) then
- deltat_int2 = 0
- else
- deltat_int2 = NINT(deltat*1.0d6, kind=2) ! deltat (unit: 10^{-6} second)
- endif
-
- ! 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)=int(NSTEP, kind=2)
- write(12,rec=(irec-1)*60+(irec-1)*NSTEP+29) header2
- header2(1)=deltat_int2
- header2(2)=0 ! dummy
- write(12,rec=(irec-1)*60+(irec-1)*NSTEP+30) header2
- 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)=int(NSTEP, kind=2)
- write(14,rec=(irec-1)*60+(irec-1)*NSTEP+29) header2
- header2(1)=deltat_int2
- header2(2)=0 ! dummy
- write(14,rec=(irec-1)*60+(irec-1)*NSTEP+30) header2
- endif
- endif
- endif
-
- ! the "60" in the following corresponds to 240 bytes header (note the reclength is 4 bytes)
- do isample = 1, seismo_current
- write(12,rec=irec*60+(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,1))
- if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then
- write(14,rec=irec*60+(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2))
- endif
- enddo
+ call write_output_SU(x_source,z_source,irec,buffer_binary,number_of_components)
endif
More information about the CIG-COMMITS
mailing list