[cig-commits] r8570 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:57:21 PST 2007
Author: walter
Date: 2007-12-07 15:57:20 -0800 (Fri, 07 Dec 2007)
New Revision: 8570
Modified:
seismo/2D/SPECFEM2D/trunk/checkgrid.F90
seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/constants.h
seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90
seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.csh
seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90
seismo/2D/SPECFEM2D/trunk/createnum_fast.f90
seismo/2D/SPECFEM2D/trunk/createnum_slow.f90
seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
- fixed a bug detected by Kasper van Wijk: top absorbing surface was deactivated by top free surface in the acoustic case
- copied the new (improved) convolution code and script from the 3D code
- fixed small bugs in the routine to check the grid when the whole mesh is acoustic: S wave information was calculated and was incorrect because there is no S wave velocity
- replaced all "end if" with "endif" and "end do" with "enddo"
Modified: seismo/2D/SPECFEM2D/trunk/checkgrid.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2007-12-07 23:57:20 UTC (rev 8570)
@@ -1318,9 +1318,16 @@
!---- compute parameters for the spectral elements
vpmin = HUGEVAL
- vsmin = HUGEVAL
vpmax = -HUGEVAL
- vsmax = -HUGEVAL
+
+ if(any_elastic) then
+ vsmin = HUGEVAL
+ vsmax = -HUGEVAL
+ else
+ vsmin = 0
+ vsmax = 0
+ endif
+
densmin = HUGEVAL
densmax = -HUGEVAL
@@ -1330,10 +1337,16 @@
courant_stability_number_max = -HUGEVAL
lambdaPmin = HUGEVAL
- lambdaSmin = HUGEVAL
lambdaPmax = -HUGEVAL
- lambdaSmax = -HUGEVAL
+ if(any_elastic) then
+ lambdaSmin = HUGEVAL
+ lambdaSmax = -HUGEVAL
+ else
+ lambdaSmin = 0
+ lambdaSmax = 0
+ endif
+
do ispec=1,nspec
material = kmato(ispec)
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2007-12-07 23:57:20 UTC (rev 8570)
@@ -297,9 +297,9 @@
if(.not. elastic(ispec_selected_source)) then
call exit_MPI('cannot have moment tensor source in acoustic element')
- end if
+ endif
endif
- end if
+ endif
else
call exit_MPI('wrong source type')
endif
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-12-07 23:57:20 UTC (rev 8570)
@@ -480,7 +480,7 @@
endif
endif
- end if
+ endif
else
call exit_MPI('wrong source type')
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-07 23:57:20 UTC (rev 8570)
@@ -78,8 +78,11 @@
! number of iterations to solve the system for xi and eta
integer, parameter :: NUM_ITER = 4
-! error function source decay rate for Heaviside
- double precision, parameter :: SOURCE_DECAY_RATE = 1.628d0
+! we mimic a triangle of half duration equal to half_duration_triangle
+! using a Gaussian having a very close shape, as explained in Figure 4.2
+! of the manual. This source decay rate to mimic an equivalent triangle
+! was found by trial and error
+ double precision, parameter :: SOURCE_DECAY_MIMIC_TRIANGLE = 1.628d0
! non linear display to enhance small amplitudes in color images
double precision, parameter :: POWER_DISPLAY_COLOR = 0.30d0
Modified: seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/construct_acoustic_surface.f90 2007-12-07 23:57:20 UTC (rev 8570)
@@ -24,7 +24,7 @@
e2 = surface(4,i)
do k = 1, ngnod
n(k) = knods(k,tab_surface(1,i))
- end do
+ enddo
call get_acoustic_edge ( ngnod, n, type, e1, e2, ixmin, ixmax, izmin, izmax )
@@ -34,7 +34,7 @@
tab_surface(5,i) = izmax
- end do
+ enddo
end subroutine construct_acoustic_surface
@@ -59,25 +59,25 @@
ixmax = 1
izmin = 1
izmax = 1
- end if
+ endif
if ( e1 == n(2) ) then
ixmin = NGLLX
ixmax = NGLLX
izmin = 1
izmax = 1
- end if
+ endif
if ( e1 == n(3) ) then
ixmin = NGLLX
ixmax = NGLLX
izmin = NGLLZ
izmax = NGLLZ
- end if
+ endif
if ( e1 == n(4) ) then
ixmin = 1
ixmax = 1
izmin = NGLLZ
izmax = NGLLZ
- end if
+ endif
else
if ( e1 == n(1) ) then
@@ -87,13 +87,13 @@
ixmax = NGLLX
izmax = 1
- end if
+ endif
if ( e2 == n(4) ) then
ixmax = 1
izmax = NGLLZ
- end if
- end if
+ endif
+ endif
if ( e1 == n(2) ) then
ixmin = NGLLX
izmin = 1
@@ -101,14 +101,14 @@
ixmax = NGLLX
izmax = NGLLZ
- end if
+ endif
if ( e2 == n(1) ) then
ixmax = ixmin
ixmin = 1
izmax = 1
- end if
- end if
+ endif
+ endif
if ( e1 == n(3) ) then
ixmin = NGLLX
izmin = NGLLZ
@@ -117,14 +117,14 @@
ixmin = 1
izmax = NGLLZ
- end if
+ endif
if ( e2 == n(2) ) then
ixmax = NGLLX
izmax = izmin
izmin = 1
- end if
- end if
+ endif
+ endif
if ( e1 == n(4) ) then
ixmin = 1
izmin = NGLLZ
@@ -133,14 +133,14 @@
izmax = izmin
izmin = 1
- end if
+ endif
if ( e2 == n(3) ) then
ixmax = NGLLX
izmax = NGLLZ
- end if
- end if
- end if
+ endif
+ endif
+ endif
end subroutine get_acoustic_edge
Modified: seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.csh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.csh 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.csh 2007-12-07 23:57:20 UTC (rev 8570)
@@ -1,19 +1,24 @@
#!/bin/csh
-set hdur = 11.2
+# we mimic a triangle of half duration equal to half_duration_triangle
+# using a Gaussian having a very close shape, as explained in Figure 4.2
+# of the manual
+set half_duration_triangle = 11.2
+
foreach file ( $* )
set nlines = `wc -l $file `
echo $nlines > input_convolve_code.txt
-echo $hdur >> input_convolve_code.txt
+echo $half_duration_triangle >> input_convolve_code.txt
# use .true. for a triangle and .false. for a Gaussian
#echo ".true." >> input_convolve_code.txt
echo ".false." >> input_convolve_code.txt
-echo convolving $file with hdur = $hdur using lines $nlines
+echo convolving $file with half_duration_triangle = $half_duration_triangle using lines $nlines
./xconvolve_source_timefunction < $file > ${file}.convolved
+
rm input_convolve_code.txt
end
Modified: seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/convolve_source_timefunction.f90 2007-12-07 23:57:20 UTC (rev 8570)
@@ -17,82 +17,102 @@
! convolve seismograms computed for a Heaviside with given source time function
!
+! we mimic a triangle of half duration equal to half_duration_triangle
+! using a Gaussian having a very close shape, as explained in Figure 4.2
+! of the manual
+
implicit none
include "constants.h"
- integer i,j,N_j
- integer number_remove
+ integer :: i,j,N_j,number_remove,nlines
+
+ double precision :: alpha,dt,tau_j,source,exponent,t1,t2,displ1,displ2,gamma,height,half_duration_triangle
+
+ logical :: triangle
+
double precision, dimension(:), allocatable :: time,sem,sem_fil
- double precision alpha,dt,tau_j,source,exponent
- double precision t1,t2,displ1,displ2,gamma,height
- integer nlines
- double precision hdur
- logical triangle
-
! read file with number of lines in input
- open(unit=33,file='input_convolve_code.txt',status='old')
+ open(unit=33,file='input_convolve_code.txt',status='old',action='read')
read(33,*) nlines
- read(33,*) hdur
+ read(33,*) half_duration_triangle
read(33,*) triangle
close(33)
! allocate arrays
allocate(time(nlines),sem(nlines),sem_fil(nlines))
- do i=1,nlines
+! read the input seismogram
+ do i = 1,nlines
read(5,*) time(i),sem(i)
enddo
- alpha=SOURCE_DECAY_RATE/hdur
- dt=time(2)-time(1)
- N_j=int(hdur/dt)
- do i=1,nlines
- sem_fil(i)=0.
- do j=-N_j,N_j
- tau_j=dble(j)*dt
+! define a Gaussian with the right exponent to mimic a triangle of equivalent half duration
+ alpha = SOURCE_DECAY_MIMIC_TRIANGLE/half_duration_triangle
+! compute the time step
+ dt = time(2) - time(1)
+
+! number of integers for which the source wavelet is different from zero
+ if(triangle) then
+ N_j = ceiling(half_duration_triangle/dt)
+ else
+ N_j = ceiling(1.5d0*half_duration_triangle/dt)
+ endif
+
+ do i = 1,nlines
+
+ sem_fil(i) = 0.d0
+
+ do j = -N_j,N_j
+
+ if(i > j .and. i-j <= nlines) then
+
+ tau_j = dble(j)*dt
+
! convolve with a triangle
if(triangle) then
- height = 1. / hdur
- if(abs(tau_j) > hdur) then
- source = 0.
- else if (tau_j < 0) then
+ height = 1.d0 / half_duration_triangle
+ if(abs(tau_j) > half_duration_triangle) then
+ source = 0.d0
+ else if (tau_j < 0.d0) then
t1 = - N_j * dt
- displ1 = 0.
- t2 = 0
+ displ1 = 0.d0
+ t2 = 0.d0
displ2 = height
gamma = (tau_j - t1) / (t2 - t1)
- source= (1. - gamma) * displ1 + gamma * displ2
+ source= (1.d0 - gamma) * displ1 + gamma * displ2
else
- t1 = 0
+ t1 = 0.d0
displ1 = height
t2 = + N_j * dt
- displ2 = 0.
+ displ2 = 0.d0
gamma = (tau_j - t1) / (t2 - t1)
- source= (1. - gamma) * displ1 + gamma * displ2
+ source= (1.d0 - gamma) * displ1 + gamma * displ2
endif
else
! convolve with a Gaussian
- exponent = alpha*alpha*tau_j*tau_j
- if(exponent < 100.) then
+ exponent = alpha**2 * tau_j**2
+ if(exponent < 50.d0) then
source = alpha*exp(-exponent)/sqrt(PI)
else
- source = 0.
+ source = 0.d0
endif
endif
- if(i > j .and. i-j <= nlines) sem_fil(i) = sem_fil(i)+sem(i-j)*source*dt
+ sem_fil(i) = sem_fil(i) + sem(i-j)*source*dt
+ endif
+
enddo
enddo
! compute number of samples to remove from end of seismograms
- number_remove = int(hdur / dt) + 1
+ number_remove = N_j + 1
do i=1,nlines - number_remove
write(*,*) sngl(time(i)),' ',sngl(sem_fil(i))
enddo
Modified: seismo/2D/SPECFEM2D/trunk/createnum_fast.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/createnum_fast.f90 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/createnum_fast.f90 2007-12-07 23:57:20 UTC (rev 8570)
@@ -198,7 +198,7 @@
! verification de la coherence de la numerotation generee
if(minval(ibool) /= 1 .or. maxval(ibool) /= npoin) then
call exit_MPI('Error while generating global numbering')
- end if
+ endif
write(IOUT,*)
write(IOUT,*) 'Total number of points of the global mesh: ',npoin
Modified: seismo/2D/SPECFEM2D/trunk/createnum_slow.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/createnum_slow.f90 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/createnum_slow.f90 2007-12-07 23:57:20 UTC (rev 8570)
@@ -241,12 +241,12 @@
! verifier que le point de depart n'existe pas deja
if(ibool(iloc,jloc,numelem) /= 0) then
call exit_MPI('point genere deux fois')
- end if
+ endif
! verifier que le point d'arrivee existe bien deja
if(ibool(i2,j2,num2) == 0) then
call exit_MPI('point inconnu dans le maillage')
- end if
+ endif
! affecter le meme numero
ibool(iloc,jloc,numelem) = ibool(i2,j2,num2)
@@ -282,7 +282,7 @@
! verification de la coherence de la numerotation generee
if(minval(ibool) /= 1 .or. maxval(ibool) /= npoin) then
call exit_MPI('Error while generating global numbering')
- end if
+ endif
write(IOUT,*) 'Total number of points of the global mesh: ',npoin
write(IOUT,*) 'distributed as follows:'
Modified: seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/enforce_acoustic_free_surface.f90 2007-12-07 23:57:20 UTC (rev 8570)
@@ -47,10 +47,10 @@
potential_dot_acoustic(iglob) = ZERO
potential_dot_dot_acoustic(iglob) = ZERO
- end do
- end do
+ enddo
+ enddo
- end do
+ enddo
end subroutine enforce_acoustic_free_surface
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-12-07 23:57:20 UTC (rev 8570)
@@ -203,14 +203,14 @@
if ( nproc <= 0 ) then
print *, 'Number of processes (nproc) must be greater than or equal to one.'
stop
- end if
+ endif
#ifndef USE_MPI
if ( nproc > 1 ) then
print *, 'Number of processes (nproc) must be equal to one when not using MPI.'
print *, 'Please recompile with -DUSE_MPI in order to enable use of MPI.'
stop
- end if
+ endif
#endif
@@ -226,7 +226,7 @@
do i = 1, 5
metis_options = iachar(partitionning_strategy(i:i)) - iachar('0')
end do
- end if
+ endif
case(3)
scotch_strategy = trim(partitionning_strategy)
@@ -247,7 +247,7 @@
print *, 'Number of control nodes must be equal to four when reading from external mesh.'
print *, 'ngnod = 9 is not yet supported.'
stop
- end if
+ endif
call read_value_logical(IIN,IGNORE_JUNK,initialfield)
call read_value_logical(IIN,IGNORE_JUNK,assign_external_model)
@@ -350,8 +350,8 @@
end do
end do
- end if
- end if
+ endif
+ endif
! read absorbing boundaries parameters
call read_value_logical(IIN,IGNORE_JUNK,any_abs)
@@ -364,7 +364,7 @@
absright = .false.
abstop = .false.
absleft = .false.
- end if
+ endif
! read time step parameters
call read_value_integer(IIN,IGNORE_JUNK,nt)
@@ -572,7 +572,7 @@
if(minval(num_material) <= 0) stop 'Velocity model not entirely set...'
- end if
+ endif
close(IIN)
@@ -716,10 +716,10 @@
end do
end do
- end if
+ endif
else
call read_nodes_coords(nodes_coords_file, nnodes, nodes_coords)
- end if
+ endif
if ( read_external_mesh ) then
@@ -728,17 +728,22 @@
if ( any_abs ) then
call read_abs_surface(absorbing_surface_file, nelemabs, abs_surface, num_start)
- end if
+ endif
else
+
! count the number of acoustic free-surface elements
nelem_acoustic_surface = 0
+
+! if the top surface is absorbing, it cannot be free at the same time
+ if(.not. abstop) then
+
j = nzread
do i = 1,nxread
imaterial_number = num_material((j-1)*nxread+i)
if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. cs(imaterial_number) < TINYVAL ) then
nelem_acoustic_surface = nelem_acoustic_surface + 1
- end if
+ endif
enddo
allocate(acoustic_surface(4,nelem_acoustic_surface))
@@ -753,8 +758,10 @@
acoustic_surface(2,nelem_acoustic_surface) = 2
acoustic_surface(3,nelem_acoustic_surface) = elmnts(3+ngnod*((j-1)*nxread+i-1))
acoustic_surface(4,nelem_acoustic_surface) = elmnts(2+ngnod*((j-1)*nxread+i-1))
- end if
+ endif
end do
+
+ endif
!
!--- definition of absorbing boundaries
@@ -779,33 +786,33 @@
abs_surface(2,nelemabs) = 2
abs_surface(3,nelemabs) = elmnts(0+ngnod*(inumelem-1))
abs_surface(4,nelemabs) = elmnts(1+ngnod*(inumelem-1))
- end if
+ endif
if(absright .and. ix == nxread) then
nelemabs = nelemabs + 1
abs_surface(1,nelemabs) = inumelem-1
abs_surface(2,nelemabs) = 2
abs_surface(3,nelemabs) = elmnts(1+ngnod*(inumelem-1))
abs_surface(4,nelemabs) = elmnts(2+ngnod*(inumelem-1))
- end if
+ endif
if(abstop .and. iz == nzread) then
nelemabs = nelemabs + 1
abs_surface(1,nelemabs) = inumelem-1
abs_surface(2,nelemabs) = 2
abs_surface(3,nelemabs) = elmnts(3+ngnod*(inumelem-1))
abs_surface(4,nelemabs) = elmnts(2+ngnod*(inumelem-1))
- end if
+ endif
if(absleft .and. ix == 1) then
nelemabs = nelemabs + 1
abs_surface(1,nelemabs) = inumelem-1
abs_surface(2,nelemabs) = 2
abs_surface(3,nelemabs) = elmnts(0+ngnod*(inumelem-1))
abs_surface(4,nelemabs) = elmnts(3+ngnod*(inumelem-1))
- end if
+ endif
end do
end do
- end if
+ endif
- end if
+ endif
! compute min and max of X and Z in the grid
@@ -876,7 +883,7 @@
print *,'Grid saved in Gnuplot format...'
print *
- end if
+ endif
!*****************************
@@ -894,14 +901,14 @@
if ( nproc > 1 ) then
call mesh2dual_ncommonnodes(nelmnts, (nxread+1)*(nzread+1), elmnts_bis, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
- end if
+ endif
else
if ( nproc > 1 ) then
call mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
- end if
+ endif
- end if
+ endif
if ( nproc == 1 ) then
@@ -941,7 +948,7 @@
end select
- end if
+ endif
! beware of fluid solid edges : coupled elements are transfered to the same partition
if ( ngnod == 9 ) then
@@ -950,7 +957,7 @@
else
call acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts, nb_materials, cs, num_material, &
nproc, part, nedges_coupled, edges_coupled)
- end if
+ endif
! local number of each element for each partition
call Construct_glob2loc_elmnts(nelmnts, part, nproc, glob2loc_elmnts)
@@ -959,7 +966,7 @@
if ( nproc > 1 ) then
deallocate(nnodes_elmnts)
deallocate(nodes_elmnts)
- end if
+ endif
allocate(nnodes_elmnts(0:nnodes-1))
allocate(nodes_elmnts(0:nsize*nnodes-1))
nnodes_elmnts(:) = 0
@@ -981,9 +988,9 @@
end do
- end if
+ endif
- end if
+ endif
! local number of each node for each partition
@@ -998,12 +1005,12 @@
else
call Construct_interfaces(nelmnts, nproc, part, elmnts, xadj, adjncy, tab_interfaces, &
tab_size_interfaces, ninterfaces, nb_materials, cs, num_material)
- end if
+ endif
print *, '04'
allocate(my_interfaces(0:ninterfaces-1))
allocate(my_nb_interfaces(0:ninterfaces-1))
print *, '05'
- end if
+ endif
! setting absorbing boundaries by elements instead of edges
if ( any_abs ) then
@@ -1014,7 +1021,7 @@
nelmnts, &
elmnts, ngnod)
print *, 'nelemabs_merge', nelemabs_merge
- end if
+ endif
! *** generate the databases for the solver
@@ -1092,7 +1099,7 @@
glob2loc_elmnts, part, iproc, 1)
else
nelemabs_loc = 0
- end if
+ endif
call Write_surface_database(15, nelem_acoustic_surface, acoustic_surface, nelem_acoustic_surface_loc, &
iproc, glob2loc_elmnts, &
@@ -1133,7 +1140,7 @@
write(15,*) 'Interfaces:'
write(15,*) 0, 0
- end if
+ endif
write(15,*) 'List of absorbing elements (bottom right top left):'
@@ -1143,7 +1150,7 @@
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right, &
glob2loc_elmnts, part, iproc, 2)
- end if
+ endif
write(15,*) 'List of acoustic free-surface elements:'
call Write_surface_database(15, nelem_acoustic_surface, acoustic_surface, nelem_acoustic_surface_loc, &
@@ -1215,7 +1222,7 @@
enddo
close(15)
- end if
+ endif
print *
@@ -1268,15 +1275,12 @@
else
num_9 = (nx+1)*(nz+1) + nx*(nz+1) + floor(real(j)/real(2))*(nx*2+1) + i + 1
- end if
- end if
+ endif
+ endif
-
end function num_9
-
-
!--- spline to describe the interfaces
double precision function value_spline(x,xinterface,zinterface,coefs_interface,npoints_interface)
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-08-29 02:20:23 UTC (rev 8569)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-07 23:57:20 UTC (rev 8570)
@@ -1360,7 +1360,7 @@
else if(time_function_type == 5) then
hdur = 1.d0 / f0
hdur_gauss = hdur * 5.d0 / 3.d0
- source_time_function(it) = factor * 0.5d0*(1.0d0+erf(SOURCE_DECAY_RATE*(time-t0)/hdur_gauss))
+ source_time_function(it) = factor * 0.5d0*(1.0d0 + erf(SOURCE_DECAY_MIMIC_TRIANGLE*(time-t0)/hdur_gauss))
else
call exit_MPI('unknown source time function')
More information about the cig-commits
mailing list