- added Brune and Smoothed Brune source time functions

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. Please note that we have implemented time derivatives of the Brune and smoothed Brune functions as the moment rate functions. For these source time functions, \texttt{hdurorf0} is the source duration or 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,85 @@ end function comp_source_time_function_d2rck
 !-------------------------------------------------------------------------------------------------
 !
 
+  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 :: omega,omegat
+
+  ! Brune source-time function
+  ! Moment function
+  !if(t .lt. 0.d0)then
+  !  comp_source_time_function_brune = 0.d0
+  !else 
+  !  omegat = 2.d0*PI*f0*t
+  !  comp_source_time_function_brune = 1.d0 - exp( -omegat ) * (1.0d0+omegat) 
+  !endif
+  ! Moment rate function
+  if(t .lt. 0.d0)then
+    comp_source_time_function_brune = 0.d0
+  else 
+    omega = 2.d0*PI*f0
+    omegat = omega*t
+    comp_source_time_function_brune = omega*omegat*exp( -omegat ) 
+  endif
+
+  end function comp_source_time_function_brune
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+  double precision function comp_source_time_function_smooth_brune(t,f0)
+
+  use constants, only: PI
+
+  implicit none
+
+  double precision, intent(in) :: t,f0
+
+  ! local variables
+  double precision,parameter :: tau0=2.31d0
+  double precision :: omega,omegat
+
+  ! Brune source-time function
+  ! Moment function
+  !omegat = 2.d0*PI*f0*t
+  !if (t .lt. 0.d0) then
+  !  comp_source_time_function_smooth_brune = 0.d0
+  !elseif (omegat .ge. 0.d0 .and. omegat .lt. tau0)then 
+  !  comp_source_time_function_smooth_brune = 1.d0 - exp(-omegat)*( 1.0d0 + omegat +  &
+  !  0.5d0*omegat**2 - (1.5d0*omegat**3)/tau0 + (1.5d0*omegat**4)/(tau0**2) - &
+  !  (0.5d0*omegat**5)/(tau0**3) )
+  !else ! (omegat .gt. tau0)then 
+  !  comp_source_time_function_smooth_brune = 1.d0 - exp( -omegat ) *
+  !  (1.0d0+omegat) 
+  !endif
+  ! Moment rate function
+  omega = 2.d0*PI*f0
+  omegat = omega*t
+  if (t .lt. 0.d0) then
+    comp_source_time_function_smooth_brune = 0.d0
+  elseif (omegat .ge. 0.d0 .and. omegat .lt. tau0)then 
+    comp_source_time_function_smooth_brune = ( 0.5d0*omega*(omegat**2)* &
+      exp(-omegat)/tau0**3 ) * ( tau0**3 - 3.d0*(tau0**2)*(omegat-3.d0) + &
+      3.d0*tau0*omegat*(omegat-4.d0) - (omegat**2)*(omegat-5.d0) )
+  else ! (omegat .gt. tau0)then 
+    comp_source_time_function_smooth_brune = omega*omegat*exp( -omegat ) 
+  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

@@ -857,6 +857,7 @@ double precision function get_stf_viscoelastic(time_source_dble,isource,it_tmp_e
 
   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 @@ double precision function get_stf_viscoelastic(time_source_dble,isource,it_tmp_e
     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 source duration or the rise time
+      ! Frequency parameter:
+      f0=1.d0/hdur(isource)
+      stf = comp_source_time_function_brune(time_source_dble,f0)
+    case (6)
+      ! Smoothed Brune source time function
+      ! hdur is the source duration or the rise time
+      ! Frequency parameter:
+      f0=1.d0/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 @@ subroutine get_force(FORCESOLUTION,tshift_src,hdur,lat,long,depth,NSOURCES, &
       ! 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 @@ subroutine locate_source()
               ! 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 @@ subroutine setup_stf_constants()
       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