Brune and Smoothed Brune source time functions

DATA/CMTSOLUTION

@@ -0,0 +1,13 @@
+PDE  1999 01 01 00 00 00.00  67000 67000 -25000 4.2 4.2 homog_test
+event name:       homog_test
+time shift:       0.0000
+half duration:    5.0
+latorUTM:       67000.0
+longorUTM:      67000.0
+depth:          30.0
+Mrr:      -7.600000e+27
+Mtt:       7.700000e+27
+Mpp:      -2.000000e+26
+Mrt:      -2.500000e+28
+Mrp:       4.000000e+26
+Mtp:      -2.500000e+27

DATA/STATIONS

@@ -0,0 +1,4 @@
+X20 DB 67000.00 22732.14 0.0 0.0
+X30 DB 67000.00 34696.43 0.0 0.0
+X40 DB 67000.00 46660.71 0.0 0.0
+X50 DB 67000.00 58625.00 0.0 0.0

doc/USER_MANUAL/05_running_the_solver.tex

@@ -193,7 +193,7 @@ \chapter{Running the Solver \texttt{xspecfem3D}}\label{cha:Running-the-Solver}
 \item \texttt{latorUTM:} Set the latitude or UTM $x$ coordinate.
 \item \texttt{longorUTM:} Set the longitude or UTM $y$ coordinate.
 \item \texttt{depth:} Set the depth of the source (in km).
-\item \texttt{source time function:} Set the type of source-time function: 0 = Gaussian function, 1 = Ricker function, 2 = Heaviside (step) function, 3 = monochromatic function, 4 = Gaussian function as defined in \citet{Meschede2011}.
+\item \texttt{source time function:} Set the type of source-time function: 0 = Gaussian function, 1 = Ricker function, 2 = Heaviside (step) function, 3 = monochromatic function, 4 = Gaussian function as defined in \citet{Meschede2011}, 5 = Brune function, and 6 = Smoothed Brune function. For the Brune and Smoothed Brune source time functions \texttt{hdurorf0} is the rise time. The Brune and Smoothed Brune functions are currently implemented only for the viscoelastic simulations.
 When {\texttt{USE\_RICKER\_TIME\_FUNCTION}} is turned on in the main parameter file \texttt{DATA/Par\_file},
 it will override this source time function type selection and always use a Ricker wavelet.
 Note that we use the standard definition of a Ricker, for a dominant frequency $f_0$:

src/specfem3D/comp_source_time_function.f90

@@ -263,6 +263,58 @@
 !-------------------------------------------------------------------------------------------------
 !
 
+  double precision function comp_source_time_function_brune(t,f0)
+
+  use constants, only: PI
+
+  implicit none
+
+  double precision, intent(in) :: t,f0
+
+  ! local variables
+  double precision :: ft
+
+  ! Brune source-time function
+  if(t .lt. 0.d0)then
+    comp_source_time_function_brune = 0.d0
+  else 
+    ft = f0*t
+    comp_source_time_function_brune = 1.d0 - exp( -ft ) * (1.0d0+ft) 
+  endif
+
+  end function comp_source_time_function_brune
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+  double precision function comp_source_time_function_smooth_brune(t,f0)
+
+  implicit none
+
+  double precision, intent(in) :: t,f0
+
+  ! local variables
+  double precision,parameter :: tau0=2.31d0
+  double precision :: f,ft
+
+  ! Brune source-time function
+  ft = f0*t
+  if (ft .lt. 0.d0) then
+    comp_source_time_function_smooth_brune = 0.d0
+  elseif (ft .ge. 0.d0 .and. ft .lt. tau0)then 
+    comp_source_time_function_smooth_brune = 1.d0 - exp(-ft)*( 1.0d0 + ft +  &
+    0.5d0*ft**2 - (1.5d0*ft**3)/tau0 + (1.5d0*ft**4)/(tau0**2) - &
+    (0.5d0*ft**5)/(tau0**3) )
+  else ! (ft .gt. tau0)then 
+    comp_source_time_function_smooth_brune = 1.d0 - exp( -ft ) * (1.0d0+ft) 
+  endif
+
+  end function comp_source_time_function_smooth_brune
+
+!
+!-------------------------------------------------------------------------------------------------
+!
 
   double precision function comp_source_time_function_mono(t,f0)
 

src/specfem3D/compute_add_sources_viscoelastic.F90

@@ -838,7 +838,7 @@
 
 ! returns source time function value for specified time
 
-  use constants, only: USE_MONOCHROMATIC_CMT_SOURCE
+  use constants, only: PI,USE_MONOCHROMATIC_CMT_SOURCE
 
   use specfem_par, only: USE_FORCE_POINT_SOURCE,USE_RICKER_TIME_FUNCTION, &
                          hdur,hdur_Gaussian,force_stf
@@ -857,6 +857,7 @@
 
   double precision, external :: comp_source_time_function,comp_source_time_function_rickr, &
     comp_source_time_function_gauss,comp_source_time_function_gauss_2, &
+    comp_source_time_function_brune,comp_source_time_function_smooth_brune, &
     comp_source_time_function_mono,comp_source_time_function_ext
 
   ! external source time function
@@ -890,6 +891,18 @@
     case (4)
       ! Gaussian source time function by Meschede et al. (2011)
       stf = comp_source_time_function_gauss_2(time_source_dble,hdur(isource))
+    case (5)
+      ! Brune source time function
+      ! hdur is the rise time
+      ! Frequency parameter:
+      f0=2.0d0*PI/hdur(isource)
+      stf = comp_source_time_function_brune(time_source_dble,f0)
+    case (6)
+      ! Smoothed Brune source time function
+      ! hdur is the rise time
+      ! Frequency parameter:
+      f0=2.0d0*PI/hdur(isource)
+      stf = comp_source_time_function_smooth_brune(time_source_dble,f0)
     case default
       stop 'unsupported force_stf value!'
     end select

src/specfem3D/get_force.f90

@@ -250,6 +250,20 @@
       ! null half-duration indicates a Dirac
       ! replace with very short Gaussian function
       if (hdur(isource) < 5. * DT ) hdur(isource) = 5. * DT
+    case (5)
+      ! Brune source time function
+      ! half-duration is the rise time
+      ! (see constants.h: TINYVAL = 1.d-9 )
+      if (hdur(isource) < TINYVAL ) then
+        stop 'Error set force period, make sure all forces have a non-zero rise time'
+      endif
+    case (6)
+      ! Smoothed Brune source time function
+      ! half-duration is the rise time
+      ! (see constants.h: TINYVAL = 1.d-9 )
+      if (hdur(isource) < TINYVAL ) then
+        stop 'Error set force period, make sure all forces have a non-zero rise time'
+      endif
     case default
       stop 'unsupported source time function type (force_stf) value!'
     end select

src/specfem3D/locate_source.F90

@@ -505,6 +505,18 @@
               ! Gaussian by Meschede et al. (2011)
               write(IMAIN,*) '    using Gaussian source time function by Meschede et al. (2011), eq.(2)'
               write(IMAIN,*) '             tau: ',hdur(isource),' seconds'
+            case (5)
+              ! Brune
+              write(IMAIN,*) '    using Brune source time function'
+              ! prints rise time for point forces
+              write(IMAIN,*)
+              write(IMAIN,*) '    rise time: ',hdur(isource),' seconds'
+            case (6)
+              ! Smoothed Brune
+              write(IMAIN,*) '    using Smoothed Brune source time function'
+              ! prints rise time for point forces
+              write(IMAIN,*)
+              write(IMAIN,*) '    rise time: ',hdur(isource),' seconds'
             case default
               stop 'unsupported force_stf value!'
             end select

src/specfem3D/setup_sources_receivers.f90

@@ -461,6 +461,14 @@
       case (4)
         ! Gaussian source time function by Meschede et al. (2011)
         t0 = min(t0,1.5d0 * (tshift_src(isource) - hdur(isource)))
+      case (5)
+        ! Brune
+        ! This needs to be CHECKED!!!
+        t0 = min(t0,1.5d0 * (tshift_src(isource) - hdur(isource)))
+      case (6)
+        ! Smotthed Brune
+        ! This needs to be CHECKED!!!
+        t0 = min(t0,1.5d0 * (tshift_src(isource) - hdur(isource)))
       case default
         stop 'unsupported force_stf value!'
       end select