[cig-commits] r8970 - seismo/2D/SPECFEM2D/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Dec 23 05:27:33 PST 2007
Author: dkomati1
Date: 2007-12-23 05:27:33 -0800 (Sun, 23 Dec 2007)
New Revision: 8970
Modified:
seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
updated several comments
Modified: seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90 2007-12-23 13:05:59 UTC (rev 8969)
+++ seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90 2007-12-23 13:27:33 UTC (rev 8970)
@@ -119,7 +119,7 @@
!!$ (t - x/2.d0 - (sqrt(3.d0)*(-9 + z))/2.d0))*exp(-(a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
!!$ (2*t - x + sqrt(3.d0)*(-9 + z))*exp(-(a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/4.d0)))/2.d0
-!to calcul the derivative of the displacement, we take the velocity ricker expression and we multiply by
+! to ompute the derivative of the displacement, we take the velocity ricker expression and we multiply by
! the derivative of the interior argument of ricker_Bielak_veloc
dxUx = A_plane(1) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) * (-sin(angleforce)/c_inc)&
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-23 13:05:59 UTC (rev 8969)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-23 13:27:33 UTC (rev 8970)
@@ -1510,21 +1510,21 @@
!=======================================================================
!
- ! Calculation of the initialfield for planewave
+ ! Calculation of the initialfield for plane wave
!
!=======================================================================
print *,'Number of grid points: ',npoin
- print *,'*** calculation of initial planewave ***'
+ print *,'*** calculation of initial plane wave ***'
if (source_type == 1) then
print *,'initial P wave of', angleforce*180.d0/pi, 'degrees introduced...'
else if (source_type == 2) then
print *,'initial SV wave of', angleforce*180.d0/pi, ' degrees introduced...'
else
- call exit_MPI('Not recognized source_type : 1 for P planewaves, 2 for SV planewaves!!!')
+ call exit_MPI('Unrecognized source_type: should be 1 for plane P waves, 2 for plane SV waves!')
endif
- !only implemented for homogeneous media so only 1 material supported
+ ! only implemented for homogeneous media therefore only 1 material supported
if (numat==1) then
mu = elastcoef(2,numat)
@@ -1534,7 +1534,7 @@
cploc = sqrt(lambdaplus2mu/denst)
csloc = sqrt(mu/denst)
- !P case
+ ! P wave case
if (source_type == 1) then
p=sin(angleforce)/cploc
@@ -1543,7 +1543,7 @@
angleforce_refl = asin(p*csloc)
- !from formulas (5.26) (5.27) p140 in Aki & Richards (1980)
+ ! from formulas (5.26) and (5.27) p 140 in Aki & Richards (1980)
PP = (- cos(2.d0*angleforce_refl)**2/csloc**3 + 4.d0*p**2*cos(angleforce)*cos(angleforce_refl)/cploc) / &
( cos(2.d0*angleforce_refl)**2/csloc**3 + 4.d0*p**2*cos(angleforce)*cos(angleforce_refl)/cploc)
@@ -1551,40 +1551,40 @@
(csloc**2*(cos(2.d0*angleforce_refl)**2/csloc**3 &
+4.d0*p**2*cos(angleforce)*cos(angleforce_refl)/cploc))
- print *,'reflected convert planewave angle: ', angleforce_refl*180.d0/pi, '\n'
+ print *,'reflected convert plane wave angle: ', angleforce_refl*180.d0/pi, '\n'
- !from Table 5.1 p141 in Aki & Richards (1980)
- !we put the oposite sign on z coefficients because z axe is oriented from bottom to top
+ ! from Table 5.1 p141 in Aki & Richards (1980)
+ ! we put the opposite sign on z coefficients because z axis is oriented from bottom to top
A_plane(1) = sin(angleforce); A_plane(2) = cos(angleforce)
B_plane(1) = PP * sin(angleforce); B_plane(2) = - PP * cos(angleforce)
C_plane(1) = PS * cos(angleforce_refl); C_plane(2) = PS * sin(angleforce_refl)
- !SV case
+ ! SV wave case
else
p=sin(angleforce)/csloc
c_inc = csloc
c_refl = cploc
- !if this coefficient is over 1, we are over the critical SV wave angle, there can't be a converted P wave
+ ! if this coefficient is greater than 1, we are beyond the critical SV wave angle and there cannot be a converted P wave
if (p*cploc<=1.d0) then
angleforce_refl = asin(p*cploc)
- !from formulas (5.30) (5.31) p140 in Aki & Richards (1980)
+ ! from formulas (5.30) and (5.31) p 140 in Aki & Richards (1980)
SS = (cos(2.d0*angleforce_refl)**2/csloc**3 - 4.d0*p**2*cos(angleforce)*cos(angleforce_refl)/cploc) / &
(cos(2.d0*angleforce_refl)**2/csloc**3 + 4.d0*p**2*cos(angleforce)*cos(angleforce_refl)/cploc)
SP = 4.d0*p*cos(angleforce)*cos(2*angleforce) / &
(cploc*csloc*(cos(2.d0*angleforce)**2/csloc**3&
+4.d0*p**2*cos(angleforce_refl)*cos(angleforce)/cploc))
- print *,'reflected convert planewave angle: ', angleforce_refl*180.d0/pi, '\n'
+ print *,'reflected convert plane wave angle: ', angleforce_refl*180.d0/pi, '\n'
else
- call exit_MPI('can t be treated for now: SV angle too high')
+ call exit_MPI('cannot be included for now: SV angle too high, beyond critical angle')
endif
- !from Table 5.1 p141 in Aki & Richards (1980)
- !we put the oposite sign on z coefficients because z axe is oriented from bottom to top
+ ! from Table 5.1 p141 in Aki & Richards (1980)
+ ! we put the opposite sign on z coefficients because z axis is oriented from bottom to top
A_plane(1) = cos(angleforce); A_plane(2) = - sin(angleforce)
B_plane(1) = SS * cos(angleforce); B_plane(2) = SS * sin(angleforce)
C_plane(1) = SP * sin(angleforce_refl); C_plane(2) = - SP * cos(angleforce_refl)
@@ -1592,7 +1592,7 @@
endif
else
- call exit_MPI('yet impossible to have several materials with planewaves')
+ call exit_MPI('not possible for now to have several materials with a plane wave (but could be done one day)')
endif
! get minimum and maximum values of mesh coordinates
@@ -1601,15 +1601,14 @@
xmax = maxval(coord(1,:))
zmax = maxval(coord(2,:))
-
- !initializing of the time offset to put the planewave not to close of the irregularity on the free surface
+ ! initialize the time offset to put the plane wave not too close to the free surface topography
if (abs(angleforce)<20.d0*pi/180.d0) then
time_offset=-1.d0*zmax/3.d0/c_inc
else
time_offset=0.d0
endif
- !to center rightly the wave
+ ! to correctly center the initial plane wave in the mesh
z0_source=zmax
x0_source=xmin + 1.d0*(xmax-xmin)/3.d0
@@ -1618,13 +1617,13 @@
x = coord(1,i)
z = coord(2,i)
- !z is from bottom to top so we take -z to make parallele with Aki & Richards
+ ! z is from bottom to top therefore we take -z to make parallel with Aki & Richards
z = z0_source - z
x = x - x0_source
t = 0.d0 + time_offset
- !formulas of the initial displacement for a planewave from Aki & Richards (1980)
+ ! formulas for the initial displacement for a plane wave from Aki & Richards (1980)
displ_elastic(1,i) = A_plane(1) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
+ B_plane(1) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ C_plane(1) * ricker_Bielak_displ(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
@@ -1632,7 +1631,7 @@
+ B_plane(2) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ C_plane(2) * ricker_Bielak_displ(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
- !formulas of the initial velocity for a planewave (first derivative in time of the displacement)
+ ! formulas for the initial velocity for a plane wave (first derivative in time of the displacement)
veloc_elastic(1,i) = A_plane(1) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
+ B_plane(1) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ C_plane(1) * ricker_Bielak_veloc(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
@@ -1640,7 +1639,7 @@
+ B_plane(2) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ C_plane(2) * ricker_Bielak_veloc(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
- !formulas of the initial acceleration for a planewave (first derivative in time of the velocity)
+ ! formulas for the initial acceleration for a plane wave (second derivative in time of the displacement)
accel_elastic(1,i) = A_plane(1) * ricker_Bielak_accel(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
+ B_plane(1) * ricker_Bielak_accel(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ C_plane(1) * ricker_Bielak_accel(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
@@ -1649,10 +1648,10 @@
+ C_plane(2) * ricker_Bielak_accel(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
enddo
- endif !add_Bielak
+ endif ! add_Bielak
write(IOUT,*) 'Max norm of initial elastic displacement = ',maxval(sqrt(displ_elastic(1,:)**2 + displ_elastic(2,:)**2))
- endif !initialfield
+ endif ! initialfield
deltatsquare = deltat * deltat
deltatcube = deltatsquare * deltat
More information about the cig-commits
mailing list