Thu Nov 6 08:29:17 PST 2014

Repository : https://github.com/geodynamics/specfem3d_globe

On branches: devel,master
Link       : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f

>---------------------------------------------------------------

commit 4063bb9a28ffebd6541398b4cf14d5943f81a97f
Author: Michael Afanasiev <mafanasiev at gmail.com>
Date:   Tue Aug 26 12:14:31 2014 +0200

    Added the option to use an external source time function. The parameter EXTERNAL_SOURCE_TIME_FUNCTION has been set in constants.h.in. If this is set to true (it is false by default), a text file called stf is looked for in the data directory. This file should contain a normalized source time function, with the correct number of timesteps, which will be multiplied by the moment tensor at each timestep. Also, if the external stf is specified, the halfduration is set to 0 (not needed), and the simulation start time is set to 0 as well.


>---------------------------------------------------------------

4063bb9a28ffebd6541398b4cf14d5943f81a97f
 DATA/Par_file                               |  2 +-
 setup/constants.h.in                        |  3 ++
 src/specfem3D/comp_source_time_function.f90 | 66 ++++++++++++++++++++++++++++-
 src/specfem3D/get_cmt.f90                   |  5 +++
 src/specfem3D/rules.mk                      |  2 +-
 src/specfem3D/setup_sources_receivers.f90   |  5 +++
 src/specfem3D/specfem3D_par.F90             |  3 ++
 7 files changed, 82 insertions(+), 4 deletions(-)

diff --git a/DATA/Par_file b/DATA/Par_file
index 789b043..b818616 100644
--- a/DATA/Par_file
+++ b/DATA/Par_file
@@ -170,7 +170,7 @@ USE_BINARY_FOR_LARGE_FILE       = .false.
 RECEIVERS_CAN_BE_BURIED         = .true.
 
 # print source time function
-PRINT_SOURCE_TIME_FUNCTION      = .false.
+PRINT_SOURCE_TIME_FUNCTION      = .true.
 
 #-----------------------------------------------------------
 #
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 2005456..187f1fd 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -262,6 +262,9 @@
 ! is located outside the mesh and therefore excluded from the station list
   double precision, parameter :: THRESHOLD_EXCLUDE_STATION = 50.d0
 
+! This parameter flags whether or not we will use an external source
+! time function. Set to false by default.
+  logical, parameter :: EXTERNAL_SOURCE_TIME_FUNCTION = .false. 
 
 !!-----------------------------------------------------------
 !!
diff --git a/src/specfem3D/comp_source_time_function.f90 b/src/specfem3D/comp_source_time_function.f90
index 45bfc76..39775a1 100644
--- a/src/specfem3D/comp_source_time_function.f90
+++ b/src/specfem3D/comp_source_time_function.f90
@@ -34,9 +34,14 @@
   double precision t,hdur
 
   double precision, external :: netlib_specfun_erf
+  double precision, external :: comp_source_time_function_external
 
-! quasi Heaviside
-  comp_source_time_function = 0.5d0*(1.0d0 + netlib_specfun_erf(t/hdur))
+  if ( EXTERNAL_SOURCE_TIME_FUNCTION ) then
+    comp_source_time_function = comp_source_time_function_external ()
+  else
+    ! quasi Heaviside
+    comp_source_time_function = 0.5d0*(1.0d0 + netlib_specfun_erf(t/hdur))
+  end if
 
   end function comp_source_time_function
 
@@ -62,3 +67,60 @@
   ! comp_source_time_function_rickr = -2.d0*PI*PI*f0*f0*f0*t * exp(-PI*PI*f0*f0*t*t)
 
   end function comp_source_time_function_rickr
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+  double precision function comp_source_time_function_external ( )
+
+  use specfem_par, only: it, stfArray_external
+  implicit none
+
+  ! On the first iteration, go get the ASCII file.
+  if ( .not. allocated (stfArray_external) ) then
+    call get_external_source_time_function ()
+  end if
+  
+  comp_source_time_function_external = stfArray_external (it)
+
+  end function comp_source_time_function_external
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+  subroutine get_external_source_time_function()
+
+  use specfem_par, only: NSTEP, stfArray_external
+  implicit none
+
+  integer :: RetCode, iterator
+  character (len=200) :: line
+
+  ! Allocate the source time function array to the number of time steps.
+  allocate ( stfArray_external (NSTEP) )
+  print *, NSTEP
+
+  ! Read in source time function.
+  open (unit=10, file='DATA/stf', status='old', form='formatted')
+
+  read_loop: do iterator=1,NSTEP
+    
+    read (10, '(A)', iostat = RetCode) line
+    
+    if ( RetCode /= 0 ) then
+      print *, "ERROR IN SOURCE TIME FUNCTION."
+      stop
+    end if
+    
+    ! Ignore lines with a hash (comments)
+    if ( index (line, "#") /= 0 ) cycle read_loop
+    
+    read (line, *) stfArray_external(iterator)
+
+  end do read_loop
+
+  close (10)
+
+  end subroutine get_external_source_time_function
diff --git a/src/specfem3D/get_cmt.f90 b/src/specfem3D/get_cmt.f90
index bffa29a..55cba78 100644
--- a/src/specfem3D/get_cmt.f90
+++ b/src/specfem3D/get_cmt.f90
@@ -145,6 +145,11 @@
 
   enddo
 
+  ! If we're using external stf, don't worry about hdur.
+  if (EXTERNAL_SOURCE_TIME_FUNCTION) then
+    hdur(:) = 0.d0
+  end if
+
   close(1)
 
   ! Sets tshift_cmt to zero to initiate the simulation!
diff --git a/src/specfem3D/rules.mk b/src/specfem3D/rules.mk
index 315cc7a..f1022bc 100644
--- a/src/specfem3D/rules.mk
+++ b/src/specfem3D/rules.mk
@@ -35,7 +35,6 @@ specfem3D_OBJECTS = \
 	$O/assemble_MPI_scalar.solver.o \
 	$O/assemble_MPI_vector.solver.o \
 	$O/comp_source_spectrum.solver.o \
-	$O/comp_source_time_function.solver.o \
 	$O/compute_adj_source_frechet.solver.o \
 	$O/convert_time.solver.o \
 	$O/define_derivation_matrices.solver.o \
@@ -51,6 +50,7 @@ specfem3D_OBJECTS = \
 # values_from_mesher.h
 specfem3D_OBJECTS += \
 	$O/asdf_data.solverstatic_module.o \
+	$O/comp_source_time_function.solverstatic.o \
 	$O/specfem3D_par.solverstatic_module.o \
 	$O/write_seismograms.solverstatic.o \
 	$O/check_stability.solverstatic.o \
diff --git a/src/specfem3D/setup_sources_receivers.f90 b/src/specfem3D/setup_sources_receivers.f90
index e1f373c..644b7d3 100644
--- a/src/specfem3D/setup_sources_receivers.f90
+++ b/src/specfem3D/setup_sources_receivers.f90
@@ -178,6 +178,11 @@
   ! define t0 as the earliest start time
   t0 = - 1.5d0*minval( tshift_cmt(:) - hdur(:) )
 
+  if ( EXTERNAL_SOURCE_TIME_FUNCTION ) then
+    hdur(:) = 0._CUSTOM_REAL
+    t0      = 0.d0
+  end if
+
   ! point force sources will start depending on the frequency given by hdur
   if( USE_FORCE_POINT_SOURCE ) then
     ! note: point force sources will give the dominant frequency in hdur,
diff --git a/src/specfem3D/specfem3D_par.F90 b/src/specfem3D/specfem3D_par.F90
index efc9893..8294c38 100644
--- a/src/specfem3D/specfem3D_par.F90
+++ b/src/specfem3D/specfem3D_par.F90
@@ -148,6 +148,9 @@ module specfem_par
   double precision, dimension(:), allocatable :: theta_source,phi_source
   double precision :: t0
 
+  ! External source time function.
+  double precision, dimension(:), allocatable :: stfArray_external
+
   !-----------------------------------------------------------------
   ! receivers
   !-----------------------------------------------------------------

