readme.txt

Frotran77编写的地震波传播在二维吸收介质的

Copyright (c) 2007 by the Society of Exploration Geophysicists.
For more information, go to http://software.seg.org/2007/2005 .
You must read and accept usage terms at:
 http://software.seg.org/disclaimer.txt  before use.

Minimal documentation on the FDTD code for vertically heterogeneous 
viscoelastic media:

Authors: Arash JafarGandomi & Hiroshi Takenaka

Date: February 2007
========================================
General code description: 

The fdtdvisco code uses a  finite difference scheme of 2nd-order accurate in time 
and 4th-order accurate in space to solve the elastic plane-wave equation for a 
vertically heterogeneous attenuative medium.
========================================
File description:

fdtdvisco.f : main file with time loop.
fdtd_vel.dat: input file for the velocity and density structure. It includes 
four columns:
	first column: depth (km)	
	second column: density (g/cm3)
	third column: S-wave velocity (km/s)
	fourth column: P-wave velocity (km/s)
The example input file contains the following values: [ 3.  2.  1.5  2.]

Note: The input velocity structure must contain both an initial wave zone 
and a computational domain zone. In the case of a downgoing wave incident 
from the upper half space, the initial wave zone is above the computational 
domain; in the case of an upgoing wave incident from the lower half space
it is beneath the computational domain.

Thickness of initial wave zone--------------------

The thickness of initial wave zone must be large enough to cover the 
total width of incident plane wave. In the case of downgoing wave it must 
be larger than or equal to the total width of incident plane wave and in the case 
of upgoing wave it must be at least equal to the total width of incident 
plane wave and absorbing zone which occupies 50 grid points.   


fdtd_relx.dat: input file for attenuation structure. It includes 3 columns:
 	first column: stress relaxation time (tausi)	
	second column: S-wave relaxation time (taues)
	third column: P-wave relaxation time (tauep)    
--------------
  There is one line row for each relaxation mechanism as follow:
         first row: stress and strain relaxation times for first relaxation 
                    mechanism             
         second row: stress and strain relaxation times for second relaxation
                    mechanism             
         third row: stress and strain relaxation times for third relaxation
                    mechanism             
         and so on ...
Note: number of relaxation mechanisms (Zener bodies) is represented by 
parameter "ln". Number of relaxations is arbitrary and it depends on the 
attenuation model.
Note: for multi-layered models the relaxation mechanisms of each layer should 
be separated by  an empty line.
Note: the initial wave zone should be homogeneous and perfectly elastic 
represented by "ln" Zener bodies (For example by putting all relaxation times 
equal to "1.0"). 

------- An example of input attenuation model for a homogeneous medium with 
Qp=Qs=50 (approximated by 5 Zener bodies)

0.256	         0.2842		     	0.2842
0.052203       0.0579			    0.0579
0.010273	     0.0114			    0.0114
0.0020218	     0.0022		  	  0.0022
0.00039789 	   0.00044166	  	0.00044166
          
fdtd_out.dat: output file of synthetic waveforms
	      first column: time (s)	
	      second column: horizontal component
	      third column: vertical component

fdtd_report.dat
This output file contains information about computation parameters.


========================================
Compilation:
The code is compiled with Fortran.
For example in Linux the following command line can be used:

f77 -o fdtdvisco.exe fdtdvisco.f

========================================

Execution:
run

./fdtdvisco.exe

to recreate the output in ./fdtd_out.dat and ./fdtd_report.dat 
========================================

Parameters: 

An example of parameters is given below with comments

fdtd_vel.dat     	may be used to input a structure model
The following parameters can be changed in the code and should then be entered 
in the code before compiling.
nt            =2000		 -->   number of time steps
dt (s)        =0.001   -->   time increment
nz            =600 		 -->   number of spatial grids
dz (km)       =0.005	 -->   spatial grid spacing
ln            =5       -->   number of relaxation mechanisms (Zener bodies)
nd            =1       -->   number of layers in the structure model including
                             lower half-space and/or upper half-space
                             (if the velocity and density values of initial 
                             wave zone are equal to those of adjacent layer
                             they can be considered as one layer)

Note: To choose an appropriate "dt" and "dz" use the relation 
dt < 0.606 dz / vmax , where vmax is the highest P-wave apparent velocity in 
the model. In the case of oblique incidence [vmax=1/(min. vertical slowness)]. 

Type of top boundary condition -------------------

itopb: 1    --> free surface
itopb: 2    --> absorbing boundary

Initial Wave position ----------------------------
The value of "source" represents the peak position of the initial wave (Ricker wavelet).It should be selected so that the whole initial wave is within the initial wave zone.
 
source:

Type of incident wave ----------------------------

ipors: 1    --> P wave
ipors: 2    --> S wave 

Propagation direction of incident waves 

ipdir:1    --> upgoing
ipdir:2    --> downgoing

Incident angle -----------------------------------
The incident angle is measured from vertical.
theta: 

Central frequency of Ricker wavelet---------------
fc:

Reference frequency of velocity structure---------
f0:

Receiver location-----------------------------------------
rcv:


==================================================
Running examples:
Keeping the default parameters as mentioned above and entering the input 
parameters as written below from the keyboard after running the code, the 
synthetic waveforms will be calculated. The results will be saved in the 
output files.
**************************************************
First Numerical Example of JafarGandomi and Takenaka (2007) in GEOPHYSICS.
Input:
the input file "fdtd_vel.dat" contains one row:" 3.  2.  1.5  2." and
"fdtd_relx.dat".
----------------   executing fdtdvisco.exe --------------------
select the type of top boundary condition (1)free surface or (2) absorbing boundary
itopb:
 2
Input the initial incident wave location 
source: 
 1.
Choose the type of incident wave:(1) for P wave or (2) for S wave
ipors: 
 1
propagation direction of incident wave (1)upgoing or (2)downgoing
ipdir:
 2
Input the incident angle (degrees measured from vertical)
theta: 
20.
Input the central frequency of Ricker wavelet
fc:
20.
Input the reference frequency of velocity structure
f0:
20
Receiver location
rcv:
 2.

Output file for waveforms:
fdtd_out.dat,    calculated waveforms for first numerical example. Includes 3 columns:
    first column: time
    second column: horizontal component                              
    third column: vertical component
                             
**************************************************
Third Numerical Example of JafarGandomi and Takenaka (2007) in GEOPHYSICS.
Input:
the input file "fdtd_vel.dat" contains one row  as "3.  2.  1.5  2.65" and
"fdtd_relx.dat".
----------------   executing fdtdvisco.exe --------------------
select the type of top boundary condition (1)free surface or (2) absorbing boundary
itopb:
 1
Input the initial incident wave location 
source: 
 2.
Choose the type of incident wave:(1) for P wave or (2) for S waves 
ipors: 
 1
propagation direction of incident waves (1)upgoing or (2)downgoing
ipdir:
 1
Input the incident angle (degrees measured from vertical)
theta: 
20.
Input the central frequency of Ricker wavelet
fc:
20.
Input the reference frequency of velocity structure
f0:
20
Receiver location 
rcv:
 0.
 
Output file for waveforms:
fdtd_out.dat,    calculated waveforms for third example. Includes 3 columns.
     first column: time
     second column: horizontal component                              
     third column: vertical component