ela2d.f
上传用户：ls4004p
上传日期：2007-08-05
资源大小：2314k
文件大小：49k
- PROGRAM ela2d
C********************************************************************
C********************************************************************
C
C EEEEEEEEE LL AAAA 2222222 DDDDDD
C EEEEEEEEE LL AAAAAA 222222222 DDDDDDD
C EE LL AAA AAA 22 222 DD DDD
C EE LL AAA AAA 222 DD DDD
C EEEEEE LL AAA AAA 222 DD DD
C EEEEEE LL AAAAAAAAAA 222 DD DD
C EE LL AAAAAAAAAA 222 DD DDD
C EE LL AA AA 222 DD DDD
C EEEEEEEEE LLLLLLLL AA AA 222222222 DDDDDDD
C EEEEEEEEE LLLLLLLL AA AA 222222222 DDDDDD
C
C********************************************************************
C
C RELEASE 1.1 (17. November 1997)
C
C********************************************************************
C Finite-Differences Modelling for Elastic Isotropic Media
C********************************************************************
C author : Joachim Falk
C Telefon (00 49)/(0) - 40 / 4123 5050
C E-mail falk@dkrz.de
C adresse : Institut fuer Geophysik , Uni Hamburg
C Bundesstrasse 55
C 20146 Hamburg / Germany
C********************************************************************
c
C FD-MODELLING for 2D hetrogeneous isotropic elastic media
C
C For the representation of acoustic media just set:
C mu = 0.0
C lamu = lambda + 2 * mu = la = lambda
C
c Features :
c
c direct solution of the elastic wave equation
c
c 6th order FD-operators for spatial derivatives
c
c 2nd order FD-time stepping formalism (leap-frog)
c or
c 4th order FD-time stepping formalism
c
c Higdon absorbing boundary condition + dissipative boundary
c
c free surface or absorbing boundary at the upper grid boundary
c
c
c !!! The program is based on a staggered grid formulation !!!
c The grid design looks like (e.g., Virieux, 1986)
c
c i-1 i-1/2 i i+1/2 i+1 i+3/2
c
c ------+---------------+----------------+----------------
c | | |
c j-1 | * + | * + | * +
c | | |
c | | |
c j-1/2 | ~ # | ~ # | ~ #
c | | |
c ------+---------------+----------------+----------------
c | | |
c j | * + | * + | * +
c | | |
c | | |
c j+1/2 | ~ # | ~ # | ~ #
c | | |
c ------+---------------+----------------+----------------
c
c Each grid cell consists of 4 points (*,+,~,#), therefore
c the grid consits of 4 grids which are staggered by half
c grid spacing. (see e.g., Virieux 1986)
c
c parameter distribution:
c
c Lame la, (la+2*mu) and stress Sxx, Szz at all '*'
c Lame mu and stress Sxz at all '#'
c density and horizontal displacement UX at all '+'
c density and vertical displacement UZ at all '~'
c
c ! user must allways consider: Ux and Uz not defined at the same
c location !
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c discretization:
c
c use at least 4-5 grid points per shortest wavelength
c (op 6th order):
c use at least 3 grid points per shortest wavelength
c (op 12th order):
c
c Lambda_min
c ---------- >= 5.0 ! or >=3. for 12th order operators
c max(dx,dz)
c
c time step dt has to be choosen sufficiently for the
c stability criterion. For minimization of numerical
c dispersion use :
c
c dt = 0.2 * dt_max for 2nd order time stepping
c
c or
c
c dt = 0.7 * dt_max for 4th order time stepping
c
c ^(-0.5)
c G { 2.48^2 2.48^2 }
c dt_max = ----- * {-------- + --------}
c c_max { dx^2 dz^2 }
c
c c_max is the highest phase velocity in the velocity model
c G = 2.0 for 2nd order time stepping
c G = 3.46 for 4th order time stepping
c CHANGE 2.48^2 to 2.68^2 for 12th order operators
c
C**************************************************************************
C
C special routines: 1) subroutines 'dxp2d_a' and 'dxm2d_a' compute
C the derivative, add the result to the OutPut field
C multiply all the Output field the MF field,
C 2) all other derivative routines just compute
C the derivative for the defined area and set Output
C field equal zero elsewhere
C
C**************************************************************************
c
c The authors appreciate to get a short reponse on
c remaining bugs as well as on useful and interesting
c features added by users.
c
C**************************************************************************
c Files:
c input files 'ela2d.in' , 'vel_P.dat' , 'vel_S.dat'
c 'density.dat' , 'grid.param'
c
c 'ela2d.in' : ASCII, setting of the input parameters
c (explanations to parameters: see source code)
C
c 'grid.param' : ASCII, contains information on nx, nz, dx, dz
C
c 'vel_P.dat', 'vel_S.dat', 'density.dat' :
c binary, each consisting of
c nx columns of nz values containing the
c velocity model [m/s] (P- and S-wave vel.)
c and density model [kg/m^3]
c
c output files 'snapX.H@' , 'snapZ,H@' , 'snapshot.H',
c 'seismicX.H@', 'seismicZ.H@' , 'seismic.H'
c 'wavelet.H@', 'wavelet.H'
c
c '*X.H@' : data (binary) (X indicates horizontal displacement)
c '*Z.H@' : data (binary) (Z indicates vertical displacement)
c '*.H' : ASCII parameter files for SEP package
C**************************************************************************
C example for
C
C01 INPUTFILE for ELA2D by Joachim Falk! DO NOT change the order of the lines
C02 Kommentar :
C03 order of time integration (2nd or 4th)___________________________________
C04 4
C05 dt[s]__________tmax[s]____________f_cent[Hz]_____________________________
C06 0.00125 0.74 25.0
C07 sp(x)________sp(z)______________________________________(source coord [m])
C08 300.0 24.0
C09 source type (1-3)_________wavelet type (1-4)_____________________________
C10 1 2
C11 active____write stride [timesteps]______________________(seismic traces)_
C12 1 2
C13 active__delay [s]_____interval [s]______________________(diashow)________
C14 1 0.1 0.08
C15 Higdon______________Free_Surf____________________________________________
C16 1 1
C17 active________taperwidth_[m]______percentage____________(tapered bound.)_
C18 1 150.0 25.0
C19 =========================================================================
C20 =========================================================================
C21 receiver-coordinates [m]
C22 10 ! number of receivers
C23 pos_X pos_Z (meter in relation to left top corner of the model)
C24 180.0000 0. Linie 1
C25 192.0000 0.
C26 204.0000 0.
C27 216.0000 0.
C28 228.0000 0.
C29 240.0000 0.
C30 252.0000 0.
C31 264.0000 0.
C32 276.0000 0.
C33 288.0000 0.
C
C
C explanations:
C lines 1-3 : comment lines
C line 4 : indicates the order of time integration to be used
C line 5 : comment line
C line 6 : discrete timestep, propagation time, center frequency
C line 7 : comment line
C line 8 : source point location
C line 9 : comment line
C line 10 : type of the source , wavelet type
C
C source type {1,2,3} 1=dilatational point source
C 2=vertically oriented single couple
C 3=horizontally oriented single couple
C
C wavelet type {1,2,3,4} 1=2nd derivative of Gauss-function
C 2=Kosloff wavelet
C 3=Ricker wavelet
C 4=Kuepper wavelet
C
C line 11 : comment line
C line 12 : jumper (seismograms on/off); write-intervall
C writestep=1 : write each sample to trace
C writestep=2 : write every 2nd sample to trace
C ...
C
C line 13 : comment line
C line 14 : jumper (snapshots on/off),
C first snapshot after "delay" seconds
C snapshot every "interval" seconds (max 10 snapshots)
C
C line 15 : comment line
C line 16 : jumper (Higdon on/off), jumper (free surface on/off),
C jumper (smoothing none/hharmonic/slowness averaging)
C if Higdon<>1 : fixed boundary (reflecting)
C if free surface<>1 : fixed or absorbing boundary
C (depends on Higdon-jumper)
C averaging: jumpsmooth=0 no smoothing
C jumpsmooth=1 harmonic aver. of the shear modulus
C jumpsmooth=2 slowness averaging
C
C line 17 : comment line
C line 18 : jumper (dissipative bound. on/off)
C width of the boundary zone [m], width should be
C approx. 2 wavelengths; source point location should
C be approx 2 wavelengths away from this zone
C percentage is a scaling factor (can be allways set
C to 25.0)
C a combination of Higdon=1 and tapeacc=1 and
C percentage=25 mosttimes leads to very accurate results
C lines 19-21: comment lines
C line 22 : number of receivers
C line 23 : comment line
C line 24-(ng+23) : x and z coordinate of ng receivers
C
C--------------------------------------------------------------------
C example for 'grid.param' :
C consisting of two lines (1 comment, 1 data)
C
C 01 -- grid.param : nx nz dx dz
C 02 120 140 5.0 5.0
C
C--------------------------------------------------------------------
C
C Compilation:
C The code is written in Fortran 77.
C param.f is an INCLUDE-file setting the dimensions of
C the arrays. If any dimension is too small, just change
C it in 'param.f' and recompile the program.
C
C use the makefile for compilation on SUN-Sparc or
C HP-Convex machines:
C
C F77COMP = f77 -c [-O3] (Sun Fortran)
C F77LOAD = f77 -o
C
C F77COMP = fc -c [-O2] (Convex Fortran)
C F77LOAD = fc -o
C
C !!! Use optimization options if available !!!
C
C**************************************************************************
C
C source code files:
C ela2d.f
C d??2d.f (routines for the spatial derivatives)
C d??2d_a.f (routines for the spatial derivatives)
C tools.f (routines for dissipative boundary)
C higdon2d.f(routines for Higdon's boundary)
C
C
C**************************************************************************
C Running the program:
C
C % ela2d.x
C
C The program reads the input files from the current
C directory. Output files are written also to the current
C directory.
C**************************************************************************
c
c References:
c
c 1) Falk, J., 1994, Diploma thesis, Univ. of Hamburg, (German language)
c
c 2) Falk, J., Tessmer, E., Gajewski, D., 1996,
c "Tube Wave Modeling by the Finite-Difference Method with
c Varying Grid Spacing", Pure and Applied Geophysics (PAGEOPH),
c in press
c
c 3) Modelling The Earth For Oil Exploration, edited by Klaus Helbig,
c 1994, Final Report of the CEC's Geoscience I Program, 447-464
c
c 4) Jastram, C., 1993, PhD thesis, Univ. of Hamburg, (German language)
c
c 5) Jastram, C., and Tessmer, E., 1994,
c "Elastic Modelling on a grid with vertically varying spacings",
c Geophysical Prospecting 42, 357-370
C
c 6) Higdon, R.L., 1991, Absorbing boundary conditions for elastic
c waves, Geophysics 65, 854-859
c
c 7) Virieux, J., 1986, P-SV wave propagation in heterogeneous
c media: a velocity stress finite-difference method,
c Geophysics 51, 889-901
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
IMPLICIT NONE
INCLUDE 'param.f'
REAL*4
& rho (MAXCOLS,MAXROWS),
& rho2 (MAXCOLS,MAXROWS),
& mu(MAXCOLS,MAXROWS), lamu(MAXCOLS,MAXROWS),
& la(MAXCOLS,MAXROWS),
& Ux(MAXCOLS,MAXROWS), Uz(MAXCOLS,MAXROWS),
& UX0(MAXCOLS,MAXROWS), UZ0(MAXCOLS,MAXROWS),
& Sxx(MAXCOLS,MAXROWS), Szz(MAXCOLS,MAXROWS),
& Sxz(MAXCOLS,MAXROWS),
& tmp3(MAXCOLS,MAXROWS),tmp4(MAXCOLS,MAXROWS)
& ,PRE(MAXCOLS,MAXROWS)
REAL*4
& wavelet(WAVESIZE), tw, pr,
& lvp(MAXROWS), rvp(MAXROWS), lvs(MAXROWS), rvs(MAXROWS),
& bvp(MAXCOLS), bvs(MAXCOLS), tvp(MAXCOLS), tvs(MAXCOLS),
& hlprow(MAXCOLS), od12, G_0, dt_max,
& dx, dz, dt, f_cent, tmax, dt2, hlpx, hlpz, rdummy,
& maxamp1, maxamp2, dx2, dz2, rmax, rmin,
& spzz, spxx, ! Quellkoordinate in Meter
& TapeCoeff(MaxTapeCoeff), ! Koeffizientenfeld getaperter Rand
& MaxDamp, ! Faktor am aeusseren Rand
& rSnapInt, rSnapDelay
C rSnapInt, rSnapDelay : first snapshot at "rSnapDelay" seconds
C snapshot interval of "rSnapInt" seconds
INTEGER
& SPX, SPZ, nx, nz, ng,
& nxm1, nxm2, nzm1, nzm2,
& nsteps, TimeInt,
& in, t, r, c, ZA, sc,
& wavelength, nSnaps, SnapInt, SnapDelay, SnapMax,
& WriteStep, TapeWidth,
& doDia, DOSeis,tapeacc, rand,
& QA, fs, nsamp, wtype, jumpsmooth
INTEGER gpos(NGMAX,2)
C_____Program-Beginning________________________________________________
OPEN(11,file='ela2d.out',form='formatted')
CLOSE(11,status='delete')
OPEN(11,file='ela2d.out',form='formatted')
WRITE(11,*)'============================================'
WRITE(11,*)' E L A 2 D . X (start)'
WRITE(11,*)'============================================'
OPEN(99,ERR=899,file='grid.param',form='formatted')
READ(99,*,ERR=899)
READ(99,*,ERR=899) nx, nz, dx, dz
CLOSE(99,ERR=899)
IF ((nx.gt.MAXCOLS) .OR. (nz.gt.MAXROWS)) THEN
WRITE(11,*)' ERROR: MAXCOLS=',maxcols
WRITE(11,*)' MAXROWS=',maxrows
WRITE(11,*)' but nx = ',nx
WRITE(11,*)' nz = ',nz
WRITE(11,*)' ! CHANGE "param.f" and recompile the program !'
GOTO 900
ENDIF
OPEN(96,file='vel_P.dat',ERR=896,
& access='direct',recl=ibyte*nz)
OPEN(97,file='vel_S.dat',ERR=897,
& access='direct',recl=ibyte*nz)
OPEN(98,file='density.dat',ERR=898,
& access='direct',recl=ibyte*nz)
DO c = 1 , nx
READ(96,ERR=886,rec=c) (lamu(c,r) , r=1, nz)
READ(97,ERR=887,rec=c) (mu(c,r) , r=1, nz)
READ(98,ERR=888,rec=c) (rho(c,r) , r=1, nz)
ENDDO
-
rmin=1.0e15
rmax=0.0
DO r = 1 , nz
DO c = 1 , nx
IF (mu(c,r) .gt. 0.05) rmin = min(rmin,mu(c,r))
rmax = max(rmax,lamu(c,r))
lamu(c,r) = lamu(c,r)*lamu(c,r)*rho(c,r)
mu(c,r) = mu(c,r)*mu(c,r)*rho(c,r)
la(c,r) = lamu(c,r)-2.*mu(c,r)
rho2(c,r) = rho(c,r)
UX(c,r) = 0.0
UZ(c,r) = 0.0
UX0(c,r) = 0.0
UZ0(c,r) = 0.0
IF ( (la(c,r).lt.0.0) .OR. (mu(c,r).lt.0.0) ) THEN
WRITE(11,*)' STOP: irregular velocity fields'
WRITE(11,*)' x = ',c,' GP'
WRITE(11,*)' z = ',r,' GP'
WRITE(11,*)' la(x,z) = ',la(c,r)
WRITE(11,*)' mu(x,z) = ',mu(c,r)
STOP
ENDIF
ENDDO
ENDDO
CLOSE(96,ERR=876)
CLOSE(97,ERR=877)
CLOSE(98,ERR=878)
WRITE(11,*)'Highest P-wave velocity in the model: ',rmax
WRITE(11,*)'Lowest S-wave velocity in the model: ',rmin
OPEN(30,ERR=830,form='formatted',file='ela2d.in')
READ(30,*,ERR=830)
READ(30,*,ERR=830)
READ(30,*,ERR=830)
READ(30,*,ERR=830) TimeInt
IF ((TimeInt.ne.2) .and. (TimeInt.ne.4)) THEN
WRITE(11,*)' Wrong parameter for the order of time '
WRITE(11,*)' integration! Set TimeInt=4th order'
TimeInt=4
ENDIF
READ(30,*,ERR=830)
READ(30,*,ERR=830) dt, tmax, f_cent
READ(30,*,ERR=830)
READ(30,*,ERR=830) spxx, spzz
READ(30,*,ERR=830)
READ(30,*,ERR=830) QA, wtype
READ(30,*,ERR=830)
READ(30,*,ERR=830) doSeis, WriteStep
READ(30,*,ERR=830)
READ(30,*,ERR=830) doDia, rsnapdelay, rsnapint
READ(30,*,ERR=830)
READ(30,*,ERR=830) rand, fs, jumpsmooth
READ(30,*,ERR=830)
READ(30,*,ERR=830) tapeacc, tw, pr
CLOSE(30,ERR=830)
-
C slowness averaging
IF(jumpsmooth.eq.2) THEN
PRINT*,' SMOOTHING : slowness averaging'
PRINT*,' and harmonic aver. of the density model'
do r=2,nz-1
do c=2,nx-1
IF ( (mu(c-1,r-1).lt.1.) .OR. (mu(c-1,r+1).lt.1.)
& .OR. (mu(c+1,r-1).lt.1.) .OR. (mu(c+1,r+1).lt.1.) )
& STOP 'ERROR mu -> 0. : averaging option not allowed!'
PRE(c,r) = 4.0 / (1./lamu(c-1,r-1) + 1./lamu(c-1,r+1)
& +1./lamu(c+1,r-1) + 1./lamu(c+1,r+1) )
rho2(c,r)= 4.0 / (1./mu(c-1,r-1) + 1./mu(c-1,r+1)
& +1./mu(c+1,r-1) + 1./mu(c+1,r+1) )
enddo
enddo
DO r=1,nz
DO c=1,nx
lamu(c,r) = PRE(c,r)*PRE(c,r)*rho(c,r)
PRE(c,r) = rho2(c,r)*rho2(c,r)*rho(c,r)
la(c,r) = lamu(c,r) - 2.0 * PRE(c,r)
PRE(c,r) = mu(c,r)
rho2(c,r) = rho(c,r)
ENDDO
ENDDO
do r=2,nz-1
do c=2,nx-1
PRE(c,r)= 4.0 / (1./mu(c ,r) + 1./mu(c ,r+1)
& +1./mu(c+1,r) + 1./mu(c+1,r+1) )
enddo
enddo
DO r=1,nz
DO c=1,nx
mu(c,r) = PRE(c,r)*PRE(c,r)*rho(c,r)
PRE(c,r) = 0.0
ENDDO
ENDDO
do r=2,nz-1
do c=2,nx-1
rho2(c,r) = 0.5 * ( rho(c,r)+rho(c,r+1) )
rho(c,r ) = 0.5 * ( rho(c,r)+rho(c+1,r) )
enddo
enddo
ENDIF
-
C_____initialize some parameters_______________________________________
nsteps = tmax / dt
IF (nsteps.gt.ntmax) THEN
WRITE(11,*) 'ERROR : number of timesteps too large!'
WRITE(11,*) ' ! Change "param.f" and recompile the program!'
GOTO 900
ENDIF
IF (nsteps/writestep.gt.tracemax) THEN
WRITE(11,*) 'ERROR : number of samples per trace too large!'
WRITE(11,*) ' ! Change "param.f" and recompile the program'
WRITE(11,*) ' or increase the write-intervall "writestep"!'
GOTO 900
ENDIF
od12 = 1./12.
dt2 = dt*dt
dx2 = dx*dx
dz2 = dz*dz
nxm1 = nx-1
nxm2 = nx-2
nzm1 = nz-1
nzm2 = nz-2
sc = 0
SnapDelay = nint(rSnapDelay*1.0 / dt)
SnapInt = nint(rSnapInt*1.0 / dt)
-
C Computing stability criterion
IF (timeint.eq.2) THEN
G_0=2.0
ELSE
G_0=3.46
ENDIF
dt_max = G_0/rmax / sqrt(2.48**2/dx2+2.48**2/dz2) ! 6th order
C dt_max = G_0/rmax / sqrt(2.68**2/dx2+2.68**2/dz2) ! 12th order
IF (dt.ge.(0.95*dt_max)) GOTO 822
IF ((timeint.eq.2) .AND. (dt.ge.(0.25*dt_max)) ) THEN
WRITE(11,*)' '
WRITE(11,*) ' WARNING : dt is larger than 0.25 * dt_max'
ENDIF
IF ((timeint.eq.4) .AND. (dt.ge.(0.75*dt_max)) ) THEN
WRITE(11,*)' '
WRITE(11,*) ' WARNING : dt is larger than 0.75 * dt_max'
ENDIF
IF ( ( (wtype.eq.1) .OR. (wtype.eq.2) )
& .AND.
& ( max(dx,dz)*5.0 .gt. (rmin/f_cent/2.0) ) )
& WRITE(11,*) ' WARNING : dx/dz may not represent
& shortest wavelength sufficiently!'
IF ( ( (wtype.eq.3) .OR. (wtype.eq.4) )
& .AND.
& ( max(dx,dz)*5.0 .gt. (rmin/f_cent/3.0) ) )
& WRITE(11,*) ' WARNING : dx/dz may not represent
& shortest wavelength sufficiently!'
WRITE(11,*)' '
C____ source coordinate
spx = INT(spxx/dx+1.5)
IF ((spx.gt.nx).or.(spx.lt.1)) THEN
WRITE(11,*)'ERROR : X-coordinate of the source point '
WRITE(11,*)' outside of the grid extensions ! '
GOTO 900
ENDIF
spz = INT(spzz/dz+1.5)
IF ((spz.gt.nz).or.(spz.lt.1)) THEN
WRITE(11,*)'ERROR : Z-coordinate of the source point '
WRITE(11,*)' outside of the grid extensions ! '
GOTO 900
ENDIF
IF (tapeacc.eq.1) THEN
tapewidth=INT(tw/dx+1.5)
IF ((tapewidth.gt.nx) .or. (tapewidth.gt.nz)) THEN
WRITE(11,*)'ERROR : width of the dissipative '
WRITE(11,*)' boundary larger than grid extensions!'
ENDIF
rdummy = tapewidth*2.*dx/
& (dt*sqrt(lamu(spx,spz)/rho(spx,spz)))
rdummy = -1.0+2.0*((pr/100.0)**(1.0/rdummy))
maxdamp=rdummy
ENDIF
ZA=2
IF (fs.eq.1) ZA=1
C_____init of the Wavelet_____________________________________________
CALL wave_5(wavelet, wavelength, f_cent, dt, wtype)
WRITE(11,*) ' '
WRITE(11,*)' wavelet length: ', wavelength, ' time steps'
IF (wavelength .gt. wavesize) THEN
WRITE(11,*) 'ERROR : status (2=error (wavelength size))'
CALL exit(2)
ENDIF
IF (QA.eq.1) THEN ! scaling of the dilatational point source
DO in = 1,wavelength
wavelet(in) = wavelet(in) * 1.0e7 *
& sqrt(2.0*la(spx,spz)+lamu(spx,spz))
ENDDO
ENDIF
open(36,access='direct',file='wavelet.H@',
& ERR=836,recl=ibyte*wavelength)
write(36,rec=1,ERR=836) (wavelet(in),in=1,wavelength)
CLOSE(36,ERR=836)
CALL wavedata(wavelength,dt,writestep)
C_____Init of the tapering zone______________________________________
IF (tapeacc.eq.1) THEN
CALL InitTapered(TapeCoeff, TapeWidth, MaxDamp)
ENDIF
C_____Init of the diashow (SnapShots)_______________________________
IF (DODia .EQ. 1) THEN
nSnaps = nint( ((tmax - rsnapdelay) / rSnapInt)) + 1
snapmax = snapdelay + (nsnaps-1)*snapint
IF (snapmax .gt. nsteps) THEN
in = mod(snapmax-1,nsteps)
in = in / snapint
snapmax = snapmax - snapint*(in + 1)
ENDIF
OPEN(21,ERR=831,FILE='snapX.H@',
& access='direct',recl=ibyte*nz)
OPEN(22,ERR=832,FILE='snapZ.H@',
& access='direct',recl=ibyte*nz)
OPEN(71,ERR=832,FILE='snapPR.H@',
& access='direct',recl=ibyte*nz)
ENDIF
C_____Init of the seismograms_______________________________________
IF (DoSeis .EQ. 1) THEN
nsamp = 0
OPEN(30,ERR=830,form='formatted',file='ela2d.in')
DO c=1,21
READ(30,*,ERR=830)
ENDDO
READ(30,*,ERR=830) ng
READ(30,*,ERR=830)
DO r = 1 , ng
READ(30,*,ERR=830) hlpx , hlpz
gpos(r,1) = INT(hlpx/dx+0.49)+1
gpos(r,2) = INT(hlpz/dz+0.49)+1
IF ( (gpos(r,1) .gt. nx) .OR.
& (gpos(r,2) .gt. nz) ) THEN
WRITE(11,*) ' ERROR : receiver location is not'
WRITE(11,*) ' within model/grid extensions!'
WRITE(11,*) ' rec No. ',r
WRITE(11,*) ' (X/Z) : ',hlpx,hlpz
GOTO 900
ENDIF
ENDDO
CLOSE(30,ERR=830)
IF (ng.eq.0) THEN
DoSeis = 0
WRITE(11,*)' NO seismograms: number of receivers = 0'
WRITE(11,*)' set DoSeis=0 !'
ELSE
OPEN(23,ERR=833,file='seismicX.H@',
& access='direct',recl=ibyte*ng)
OPEN(24,ERR=834,file='seismicZ.H@',
& access='direct',recl=ibyte*ng)
OPEN(73,ERR=834,file='seismicPR.H@',
& access='direct',recl=ibyte*ng)
nsamp = nsamp + 1
write(23,ERR=823,rec=nsamp)
& (Ux(gpos(in,1),gpos(in,2)), in=1,ng)
write(24,ERR=824,rec=nsamp)
& (Uz(gpos(in,1),gpos(in,2)), in=1,ng)
write(73,ERR=824,rec=nsamp)
& (Uz(gpos(in,1),gpos(in,2)), in=1,ng)
ENDIF
ENDIF
C_____determination of the coefficients for the absorbing boundary
C (Higdon)
IF (rand.eq.1)
& CALL InitBound(lamu,mu,rho,dx,dz,
& tvp,tvs,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz,dt,fs)
C harmonic averaging of the shear modulus
IF (jumpsmooth.eq.1) THEN
PRINT*,' SMOOTHING : harmonic averaging of the shear modulus'
PRINT*,' and density model '
do r=1,nz-1
do c=1,nx-1
IF ( (mu(c ,r).lt.1.) .OR. (mu(c ,r+1).lt.1.)
& .OR. (mu(c+1,r).lt.1.) .OR. (mu(c+1,r+1).lt.1.) )
& STOP 'ERROR mu -> 0. : averaging option not allowed !'
mu(c,r) = 4.0 / ( 1./mu(c,r) +1./mu(c,r+1)
& + 1./mu(c+1,r)+1./mu(c+1,r+1))
rho2(c,r) = 0.5 * ( rho(c,r)+rho(c,r+1) )
rho(c,r) = 0.5 * ( rho(c,r)+rho(c+1,r) )
enddo
enddo
ENDIF
-
C scale density (for iteration loop)
do r=1,nz
do c=1,nx
rho(c,r) = dt2 / rho(c,r)
rho2(c,r) = dt2 / rho2(c,r)
enddo
enddo
C_____print some parameters to standard output________________________
WRITE(11,*) 'grid rows ', nz
WRITE(11,*) 'grid columns ', nx
WRITE(11,*) 'source point [GP] ', spx , spz
WRITE(11,*) ' [m] ', spxx,spzz
WRITE(11,*) 'number of receivers ', ng
WRITE(11,*) 'source type ', QA
WRITE(11,*) ' ( 1=dilatational source'
WRITE(11,*) ' 2=vertical displacement'
WRITE(11,*) ' 3=horizontal displacement )'
WRITE(11,*) 'wavelet type ', wtype
WRITE(11,*) ' ( 1=2nd deriv. of Gauss-function'
WRITE(11,*) ' 2=Kosloff wavelet'
WRITE(11,*) ' 3=Ricker wavelet'
WRITE(11,*) ' 4=Kuepper wavelet )'
WRITE(11,*) 'free surface (1=on) ', fs
WRITE(11,*) 'first row of comput. ', ZA
WRITE(11,*) 'dx between cols 1&2 ', dx
WRITE(11,*) 'dz between lines 1&2 ', dz
WRITE(11,*) 'central frequency ', f_cent
IF (tapeacc.eq.1) THEN
WRITE(11,*) 'dissipative boundary-zone ', TapeWidth,' [GP]'
WRITE(11,*) ' or ', tw,' [m]'
WRITE(11,*) 'max. taper factor ', MaxDamp
ENDIF
WRITE(11,*) 'Higdon-Boundary-Cond.' ,rand,' 1=on'
WRITE(11,*) 'write SeisSample each', WriteStep,' Timestep'
WRITE(11,*) 'ZeitIntegration of ', TimeInt,'th order'
WRITE(11,*) 'number of time steps ', nsteps
WRITE(11,*) 'runtime ', tmax,' s'
WRITE(11,*) 'time-intervall dt ', dt,' s'
WRITE(11,*) ' '
WRITE(11,*)' SnapShots : '
WRITE(11,*)' nSnaps = ',nsnaps ,' SnapShots '
WRITE(11,*)' Delay = ',snapdelay,' timesteps'
WRITE(11,*)' or ',rsnapdelay,' seconds'
WRITE(11,*)' Intervall = ',SnapInt ,' timesteps'
WRITE(11,*)' or ',rsnapint,' seconds'
WRITE(11,*)' snapmax = ',Snapmax ,' timesteps'
C_____BEGIN time loop iteration____________________________________
WRITE(11,*) ' '
WRITE(11,*) ' !!! BEGIN of the time loop !!!'
WRITE(11,*) ' '
CALL FLUSH(11)
-
DO t=1,nsteps
IF (mod(t,25) .EQ. 0) THEN
WRITE(11,*) 'TimeStep ', t, ' of ', nsteps
CALL FLUSH(11)
ENDIF
C_____single force (vertical)__________________________________________
IF ((QA.eq.2) .AND. (t.le.wavelength)) THEN
Uz(spx,spz) = Uz(spx,spz) + wavelet(t)
ENDIF
C_____single force (horizontal)________________________________________
IF ((QA.eq.3) .AND. (t.le.wavelength)) THEN
Ux(spx,spz) = Ux(spx,spz) + wavelet(t)
ENDIF
C_____computing strain tensor__________________________________________
CALL DZm2D(Uz,Szz,nx,nz,dz,ZA,nzm1,2,nxm1,fs)
CALL DXm2D(Ux,Sxx,nx,nz,dx,ZA,nzm1,2,nxm1)
IF (fs.eq.1) THEN
DO in = 1, nx
Szz(in,1) = -la(in,1)/lamu(in,1) * Sxx(in,1)
ENDDO
ENDIF
CALL DZp2D( Ux,Sxz, nx,nz,dz,1,nzm1,1,nxm1,fs)
CALL DXp2d_a(Uz,Sxz,mu,nx,nz,dx,1,nzm1,1,nxm1)
C_____computing stress tensor__________________________________________
DO r = ZA,nzm1
DO c=2,nxm1
hlprow(c) = Sxx(c,r)
ENDDO
DO c=2,nxm1
Sxx(c,r) = lamu(c,r)*hlprow(c)+ la(c,r) * Szz(c,r)
Szz(c,r) = lamu(c,r)*Szz(c,r) + la(c,r) * hlprow(c)
ENDDO
ENDDO
IF ( (DoDia .EQ. 1) .AND.
& (t.ge.SnapDelay) .AND. (t.le.SnapMax) .AND.
& (mod(t-SnapDelay,SnapInt) .eq. 0)) THEN
DO c=1,nx
WRITE(71,ERR=821,rec=sc*nx+c)
& ((Sxx(c,r)+Szz(c,r))/(lamu(c,r)/rho(c,r)),r=1,nz)
ENDDO
ENDIF
IF ( (DoSeis .EQ. 1) .AND.
& (MOD(t, WriteStep) .EQ. 0)) THEN
nsamp = nsamp + 1
write(73,ERR=823,rec=nsamp)
& ( (Sxx(gpos(in,1),gpos(in,2)) +
& Szz(gpos(in,1),gpos(in,2)))
& /(lamu(gpos(in,1),gpos(in,2))/
& rho(gpos(in,1),gpos(in,2)) ), in=1,ng)
ENDIF
C_____explicit free surface condition_________________________________
IF (fs.eq.1) THEN
DO c = 1, nx
Szz(c,1) = 0.0
ENDDO
ENDIF
C_____dilatational source_____________________________________________
IF ((QA.eq.1) .AND. (t .LE. wavelength)) THEN
Sxx(spx,spz) = Sxx(spx,spz) + wavelet(t)
Szz(spx,spz) = Szz(spx,spz) + wavelet(t)
ENDIF
CALL DZm2D( Sxz,tmp4, nx,nz,dz,ZA,nzm1,2,nxm2,fs)
CALL DXp2d_a(Sxx,tmp4,rho, nx,nz,dx,ZA,nzm1,2,nxm2)
CALL DZp2D( Szz,tmp3, nx,nz,dz,ZA,nzm2,2,nxm1,fs)
CALL DXm2D_A(Sxz,tmp3,rho2,nx,nz,dx,ZA,nzm2,2,nxm1)
C_____Higdon_UX_1_____________________________________________________
IF (rand.eq.1) THEN
IF (fs.eq.1) THEN
CALL higfsX1(Ux,UX0,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
ELSE
CALL higX1(Ux,UX0,tvp,tvs,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
ENDIF
ENDIF
C_____Higdon_UZ_1_____________________________________________________
IF (rand.eq.1) THEN
IF (fs.eq.1) THEN
CALL higfsZ1(Uz,UZ0,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
ELSE
CALL higZ1(Uz,UZ0,tvp,tvs,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
ENDIF
ENDIF
C_____second order integration (leap frog scheme)_____________________
IF (timeint.eq.2) THEN
DO r=ZA,nzm2
DO c=2,nxm2
hlpx = Ux(c,r)
Ux(c,r) = 2.0*Ux(c,r) - UX0(c,r) +
& tmp4(c,r)
UX0(c,r) = hlpx
hlpz = Uz(c,r)
Uz(c,r) = 2.0*Uz(c,r) - UZ0(c,r) +
& tmp3(c,r)
UZ0(c,r) = hlpz
ENDDO
hlpz = Uz(nxm1,r)
Uz(nxm1,r) = 2.0*Uz(nxm1,r) - UZ0(nxm1,r) +
& tmp3(nxm1,r)
UZ0(nxm1,r) = hlpz
ENDDO
r = nzm1
DO c=2,nxm2
hlpx = UX(c,r)
Ux(c,r) = 2.0*Ux(c,r) - UX0(c,r) +
& tmp4(c,r)
UX0(c,r) = hlpx
ENDDO
ELSE
C_____(corresponds to forth order integration)
DO r=ZA,nzm2
DO c=2,nxm2
UX0(c,r) = 2.0*Ux(c,r) - UX0(c,r) +
& tmp4(c,r)
UZ0(c,r) = 2.0*Uz(c,r) - UZ0(c,r) +
& tmp3(c,r)
ENDDO
UZ0(nxm1,r) = 2.0*Uz(nxm1,r) - UZ0(nxm1,r) +
& tmp3(nxm1,r)
ENDDO
r = nzm1
DO c=2,nxm2
UX0(c,r) = 2.0*Ux(c,r) - UX0(c,r) +
& tmp4(c,r)
ENDDO
C_____computing spatial derivatives second time (4th order time integ.)
CALL DZm2D(tmp3,Szz,nx,nz,dz,ZA,nzm1,2,nxm1,fs)
CALL DXm2D(tmp4,Sxx,nx,nz,dx,ZA,nzm1,2,nxm1)
IF (fs.eq.1) THEN
DO c = 1, nx
Szz(c,1) = -la(c,1) / lamu(c,1) * Sxx(c,1)
ENDDO
ENDIF
CALL DZp2D( tmp4,Sxz, nx,nz,dz,1,nzm1,1,nxm1,fs)
CALL DXp2d_a(tmp3,Sxz,mu,nx,nz,dx,1,nzm1,1,nxm1)
DO r=1,nzm1
DO c=1,nxm1
hlprow(c) = Sxx(c,r)
ENDDO
DO c=2,nxm1
Sxx(c,r) = lamu(c,r) * hlprow(c) + la(c,r) * Szz(c,r)
Szz(c,r) = lamu(c,r) * Szz(c,r) + la(c,r) * hlprow(c)
ENDDO
ENDDO
C_____explicit free surface condition____________________________
IF (fs.eq.1) THEN
DO c = 1, nx
Szz(c,1) = 0.0
ENDDO
ENDIF
CALL DZm2D( Sxz,tmp4, nx,nz,dz,ZA,nzm1,2,nxm2,fs)
CALL Dxp2d_a(Sxx,tmp4,rho, nx,nz,dx,ZA,nzm1,2,nxm2)
CALL DZp2D( Szz,tmp3, nx,nz,dz,ZA,nzm2,2,nxm1,fs)
CALL DXm2D_A(Sxz,tmp3,rho2,nx,nz,dx,ZA,nzm2,2,nxm1)
DO r=ZA,nzm2
DO c=2,nxm2
hlpx = Ux(c,r)
Ux(c,r) = UX0(c,r) + od12 * tmp4(c,r)
UX0(c,r) = hlpx
hlpz = Uz(c,r)
Uz(c,r) = UZ0(c,r) + od12 * tmp3(c,r)
UZ0(c,r) = hlpz
ENDDO
hlpz = Uz(nxm1,r)
Uz(nxm1,r) = UZ0(nxm1,r)+od12*tmp3(nxm1,r)
UZ0(nxm1,r) = hlpz
ENDDO
r = nzm1
DO c=2,nxm2
hlpx = Ux(c,r)
Ux(c,r) = UX0(c,r) + od12 * tmp4(c,r)
UX0(c,r) = hlpx
ENDDO
ENDIF ! IF 2nd order or 4th order time integration
C_____Higdon_UX_2____________________________________________________
IF (rand.eq.1) THEN
IF (fs.eq.1) THEN
CALL higfsX2(Ux,UX0,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
ELSE
CALL higX2(Ux,UX0,tvp,tvs,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
ENDIF
ENDIF
C_____Higdon_UZ_2_____________________________________________________
IF (rand.eq.1) THEN
IF (fs.eq.1) THEN
CALL higfsZ2(Uz,UZ0,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
ELSE
CALL higZ2(Uz,UZ0,tvp,tvs,bvp,bvs,lvp,lvs,rvp,rvs,nx,nz)
ENDIF
ENDIF
IF (fs.ne.1) THEN
Ux(1,1) = 0.
Ux(nxm1,1) = 0.
Ux(nx,1) = 0.
Uz(1,1) = 0.
Uz(nx,1) = 0.
ENDIF
C_____disspative boundary____________________________________________
IF (tapeacc.eq.1) THEN
IF (fs.eq.1) THEN
CALL Taper(Ux , Uz , TapeCoeff, nx, nz, TapeWidth)
CALL Taper(UX0, UZ0, TapeCoeff, nx, nz, TapeWidth)
ELSE
CALL Taper42d(Ux , Uz ,TapeCoeff,nx, nz, TapeWidth)
CALL Taper42d(UX0, UZ0,TapeCoeff,nx, nz, TapeWidth)
ENDIF
ENDIF
C_____saving snapshots________________________________________
IF ( (DoDia .EQ. 1) .AND.
& (t.ge.SnapDelay) .AND. (t.le.SnapMax) .AND.
& (mod(t-SnapDelay,SnapInt) .eq. 0)) THEN
maxamp1 = -1.
maxamp2 = -1.
DO r=1, nz
DO c=1,nx
maxamp1 = MAX(abs(Ux(c,r)), maxamp1)
maxamp2 = MAX(abs(Uz(c,r)), maxamp2)
ENDDO
ENDDO
DO c=1,nx
WRITE(21,ERR=821,rec=sc*nx+c) (Ux(c,r),r=1,nz)
WRITE(22,ERR=822,rec=sc*nx+c) (Uz(c,r),r=1,nz)
ENDDO
sc = sc + 1
WRITE(11,*) ' '
WRITE(11,*) ' SNAPSHOT No.',sc,' : max. amplitudes (Ux, Uz)'
WRITE(11,*) ' ',MaxAmp1,MaxAmp2
WRITE(11,*) ' at timestep ', t,' of ',nsteps
WRITE(11,*) ' '
ENDIF
C_____saving seismograms______________________________________
IF ( (DoSeis .EQ. 1) .AND.
& (MOD(t, WriteStep) .EQ. 0)) THEN
write(23,ERR=823,rec=nsamp)
& (Ux(gpos(in,1),gpos(in,2)), in=1,ng)
write(24,ERR=824,rec=nsamp)
& (Uz(gpos(in,1),gpos(in,2)), in=1,ng)
ENDIF
ENDDO !!!!! END time loop !!!!!
WRITE(11,*) ' '
WRITE(11,*) ' !!! END of the time loop !!!'
WRITE(11,*) ' '
IF (DoDia .EQ. 1) THEN
CLOSE(21)
CLOSE(22)
CLOSE(71)
WRITE(11,*)' CLOSED snapshot files !'
ENDIF
IF (DoSeis.EQ. 1) THEN
CLOSE(23)
CLOSE(24)
CLOSE(73)
WRITE(11,*)' CLOSED seismic files !'
ENDIF
IF ((DoSeis.EQ.1) .AND. (nsamp.gt.0))
& CALL demux(nsamp,ng,dt,writestep)
WRITE(11,*) nsamp ,' samples have been written to seismic traces'
IF (DoDia.eq.1)
& CALL snapdata(nx,nz,dx,dz,nsnaps,rsnapdelay,rsnapint)
WRITE(11,*)'============================================'
WRITE(11,*)' E L A 2 D . X (END)'
WRITE(11,*)'============================================'
GOTO 900
C==================================================================
821 WRITE(11,*)' ERROR 821 : Write on file No. 21 (snapX.H@)!'
GOTO 900
822 WRITE(11,*)' ERROR 822 : Write on file No. 22 (snapZ.H@)!'
GOTO 900
823 WRITE(11,*)' ERROR 823 : Write on file No. 23 (seismicX.H@)!'
GOTO 900
824 WRITE(11,*)' ERROR 824 : Write on file No. 24 (seismicZ.H@)!'
GOTO 900
830 WRITE(11,*)' ERROR 830 : operation (OPEN, READ or CLOSE) '
WRITE(11,*)' on file No. 30 (ela2d.in) was not succesful'
GOTO 900
831 WRITE(11,*)' ERROR 831 : OPEN file No. 21 (snapX.H@)!'
GOTO 900
832 WRITE(11,*)' ERROR 832 : OPEN file No. 22 (snapZ.H@)!'
GOTO 900
833 WRITE(11,*)' ERROR 833 : OPEN file No. 23 (seismicX.H@)!'
GOTO 900
834 WRITE(11,*)' ERROR 834 : OPEN file No. 24 (seismicZ.H@)!'
GOTO 900
836 WRITE(11,*)' ERROR 836 : Write on file No. 36 (wavelet.H@)!'
GOTO 900
876 WRITE(11,*)' ERROR 876 : operation (CLOSE)'
WRITE(11,*)' on file No. 96 (vel_P.dat) was not succesful'
GOTO 900
877 WRITE(11,*)' ERROR 877 : operation (CLOSE)'
WRITE(11,*)' on file No. 97 (vel_S.dat) was not succesful'
GOTO 900
878 WRITE(11,*)' ERROR 878 : operation (CLOSE)'
WRITE(11,*)' on file No. 98 (density.dat) was not succesful'
GOTO 900
886 WRITE(11,*)' ERROR 886 : operation (READ)'
WRITE(11,*)' on file No. 96 (vel_P.dat) was not succesful'
GOTO 900
887 WRITE(11,*)' ERROR 887 : operation (READ)'
WRITE(11,*)' on file No. 97 (vel_S.dat) was not succesful'
GOTO 900
888 WRITE(11,*)' ERROR 888 : operation (READ)'
WRITE(11,*)' on file No. 98 (density.dat) was not succesful'
GOTO 900
896 WRITE(11,*)' ERROR 896 : operation (OPEN)'
WRITE(11,*)' on file No. 96 (vel_P.dat) was not succesful'
GOTO 900
897 WRITE(11,*)' ERROR 897 : operation (OPEN)'
WRITE(11,*)' on file No. 97 (vel_S.dat) was not succesful'
GOTO 900
898 WRITE(11,*)' ERROR 898 : operation (OPEN)'
WRITE(11,*)' on file No. 98 (density.dat) was not succesful'
GOTO 900
899 WRITE(11,*)' ERROR 899 : operation (OPEN, READ or CLOSE)'
WRITE(11,*)' on file No. 99 (grid.param) was not succesful'
GOTO 900
900 STOP
END
C____&___1_________2_________3_________4_________5_________6_________7__
SUBROUTINE WAVE_5(WAVELT,NWAVLT,F_C,DT,kind)
INCLUDE 'param.f'
REAL*4 WAVELT(wavesize)
PARAMETER ( N=2, IBEG=1)
PI = 4. * ATAN (1.)
PI2=2.*PI
IF (kind.eq.1) THEN
C WRITE(11,*)' wavelet type: Gauss'
FP = F_C*2.0/3.6
DELAY = .8/(FP)
A1 = PI2*FP*EXP(0.5)
A2 = -0.5*(PI2*FP)**2
nwavlt=wavesize
DO nn=1,Nwavlt
TAU = (NN-1)*DT-DELAY
IF(TAU.GT.DELAY) RETURN
TEMP = A2*TAU*TAU
AMPL = A1*2.
WAVELT(NN) = AMPL*(0.5 + TEMP )*EXP(TEMP)
ENDDO
GOTO 901
ENDIF
IF (kind.eq.2) THEN
C WRITE(11,*)' wavelet type: Kosloff'
AGAUSS=F_C
TCUT=1.5/AGAUSS
NWAVLT=INT(TCUT*2./DT)
NW=NWAVLT/2
DO I=1,NWAVLT
S=(-NW+I-1)*DT*AGAUSS
WAVELT(I)=EXP(-2.*S*S)*COS(PI2*S)
ENDDO
GOTO 901
ENDIF
IF (kind.eq.3) THEN
C WRITE(11,*)' wavelet type: Ricker'
sPId2 = sqrt(PI) / 2.0
cf = f_c
tp = 1.0 / f_c
ts = 1.0 / cf
NWAVLT=int(TP*2 / DT)
a = (PI * (DT-ts) / tp)
a = a*a
wavelt(1) = sPId2 * (a-0.5) * exp(-a)
DO i=2,nwavlt
a = (PI * (i*DT-ts) / tp)
a = a*a
wavelt(i) = sPId2 * (a-0.5) * exp(-a)
ENDDO
GOTO 901
ENDIF
IF (kind.eq.4) THEN
C WRITE(11,*)' wavelet type: Kuepper'
T = 1.0 / f_c
FM = ( FLOAT(N) + 2. ) / FLOAT(N)
D = N * PI / T
IEL = INT ( T/DT ) + 1
IF ( IEL.GT.wavesize ) THEN
WRITE(11,*) ' wavelet has too many samples
& ! EXIT aus Kuepper '
STOP 'Kuepper: wavelet has too many samples!'
ENDIF
DO I = 1,IEL
Wavelt(IBEG+I-1) = SIN(D*FLOAT(I-1)*DT) -
& 1./FM*SIN(FM*D*FLOAT(I-1)*DT)
ENDDO
nwavlt = IEL
GOTO 901
ELSE
WRITE(11,*)'wavelet type No. ',kind,
& ' : wrong identifier for wavelet type ! '
STOP 'wrong identifier for wavelet type ! '
ENDIF
901 RETURN
END
C================================================================
SUBROUTINE wavedata(ns, dt, ws)
IMPLICIT NONE
INCLUDE 'param.f'
INTEGER ns, ws
REAL*4 dt, dd
CHARACTER*40 text
OPEN(10,form='formatted',file='wavelet.H')
text=' in="wavelet.H@"'
WRITE(10,*) text
text=' title="ELA2D wavelet"'
WRITE(10,*) text
IF (ns.lt.10) THEN
WRITE(10,'(A4,I1)') ' n1=',ns
ELSEIF (ns.lt.100) THEN
WRITE(10,'(A4,I2)') ' n1=',ns
ELSEIF (ns.lt.1000) THEN
WRITE(10,'(A4,I3)') ' n1=',ns
ELSEIF (ns.lt.10000) THEN
WRITE(10,'(A4,I4)') ' n1=',ns
ELSEIF (ns.lt.100000) THEN
WRITE(10,'(A4,I5)') ' n1=',ns
ENDIF
WRITE(10,'(A4,I1)') ' n2=',1
WRITE(10,'(A4,I1)') ' n3=',1
WRITE(10,*) ' esize=4'
WRITE(10,*) ' o1=0'
WRITE(10,*) ' o2=1'
WRITE(10,*) ' o3=1'
WRITE(10,*) ' d3=1'
WRITE(10,*) ' d2=1'
dd = dt
IF (dd.lt.1.0) THEN
WRITE(10,'(A4,F9.8)') ' d1=',dd
ELSEIF (dd.lt.10.0) THEN
WRITE(10,'(A4,F9.7)') ' d1=',dd
ELSEIF (dd.lt.100.0) THEN
WRITE(10,'(A4,F9.6)') ' d1=',dd
ENDIF
WRITE(10,*) ' label1="sec"'
WRITE(10,*) ' label2=" "'
WRITE(10,*) ' label3=" "'
WRITE(10,*) ' data_format="xdr_float"'
CLOSE(10)
RETURN
END
C____&___1_________2_________3_________4_________5_________6_________7__
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
SUBROUTINE demux(nsamp, ng, dt, ws)
IMPLICIT NONE
INCLUDE 'param.f'
INTEGER i,j,nsamp,ng,ws
REAL*4 sdata(tracemax,ngmax), dt
CHARACTER*40 text
IF (nsamp.gt.ntmax) THEN
WRITE(11,*)'The number of smaples per trace is larger'
WRITE(11,*)' than ', ntmax,' ! Increase the increment'
WRITE(11,*)' "write each step" in "ela2d.in" or increase'
WRITE(11,*)' the parameter "ntmax" in "param.f" and'
WRITE(11,*)' recompile the source-code!'
STOP
ENDIF
C seis UX_________________________________________________________
OPEN(23,access='direct', file='seismicX.H@',
& recl=ibyte*ng)
DO i=1,nsamp
READ(23,rec=i) (sdata(i,j) , j=1,ng)
ENDDO
CLOSE(23,status='delete')
OPEN(23,access='direct', file='seismicX.H@',
& recl=ibyte*nsamp)
DO i=1,ng
DO j=2,nsamp
sdata(j,i) = sdata(j,i) + sdata(j-1,i)
ENDDO
C DO j=2,nsamp
C sdata(j,i) = sdata(j,i) + sdata(j-1,i)
C ENDDO
WRITE(23,rec=i) (sdata(j,i) , j=1,nsamp)
ENDDO
CLOSE(23)
C seis UZ_________________________________________________________
OPEN(24,access='direct', file='seismicZ.H@',
& recl=ibyte*ng)
DO i=1,nsamp
READ(24,rec=i) (sdata(i,j) , j=1,ng)
ENDDO
CLOSE(24,status='delete')
OPEN(24,access='direct', file='seismicZ.H@',
& recl=ibyte*nsamp)
DO i=1,ng
DO j=2,nsamp
sdata(j,i) = sdata(j,i) + sdata(j-1,i)
ENDDO
C DO j=2,nsamp
C sdata(j,i) = sdata(j,i) + sdata(j-1,i)
C ENDDO
WRITE(24,rec=i) (sdata(j,i) , j=1,nsamp)
ENDDO
CLOSE(24)
C seis PR_________________________________________________________
OPEN(73,access='direct', file='seismicPR.H@',
& recl=ibyte*ng)
DO i=1,nsamp
READ(73,rec=i) (sdata(i,j) , j=1,ng)
ENDDO
CLOSE(73,status='delete')
OPEN(73,access='direct', file='seismicPR.H@',
& recl=ibyte*nsamp)
DO i=1,ng
DO j=2,nsamp
sdata(j,i) = sdata(j,i) + sdata(j-1,i)
ENDDO
C DO j=2,nsamp
C sdata(j,i) = sdata(j,i) + sdata(j-1,i)
C ENDDO
WRITE(73,rec=i) (sdata(j,i) , j=1,nsamp)
ENDDO
CLOSE(73)
OPEN(10,form='formatted',file='seismic.H')
text=' in="seismic.H@"'
WRITE(10,*) text
text=' title="ELA2D seismic data"'
WRITE(10,*) text
IF (nsamp.lt.10) THEN
WRITE(10,'(A4,I1)') ' n1=',nsamp
ELSEIF (nsamp.lt.100) THEN
WRITE(10,'(A4,I2)') ' n1=',nsamp
ELSEIF (nsamp.lt.1000) THEN
WRITE(10,'(A4,I3)') ' n1=',nsamp
ELSEIF (nsamp.lt.10000) THEN
WRITE(10,'(A4,I4)') ' n1=',nsamp
ENDIF
IF (ng.lt.10) THEN
WRITE(10,'(A4,I1)') ' n2=',ng
ELSEIF (ng.lt.100) THEN
WRITE(10,'(A4,I2)') ' n2=',ng
ELSEIF (ng.lt.1000) THEN
WRITE(10,'(A4,I3)') ' n2=',ng
ELSEIF (ng.lt.10000) THEN
WRITE(10,'(A4,I4)') ' n2=',ng
ENDIF
WRITE(10,*) ' n3=1'
WRITE(10,*) ' esize=4'
WRITE(10,*) ' o1=0'
WRITE(10,*) ' o2=1'
WRITE(10,*) ' o3=1'
WRITE(10,'(A4,F10.9)') ' d1=0',(dt*ws)
WRITE(10,*) ' d2=1'
WRITE(10,*) ' d3=1'
WRITE(10,*) ' label1="sec"'
WRITE(10,*) ' label2="trace No."'
WRITE(10,*) ' label3="no label defined"'
WRITE(10,*) ' data_format="xdr_float"'
CLOSE(10)
RETURN
END
C____&___1_________2_________3_________4_________5_________6_________7__
SUBROUTINE snapdata(nx, nz, dx, dz, nsnap, rsnap, dsnap)
IMPLICIT NONE
INCLUDE 'param.f'
INTEGER nx, nz, nsnap
REAL*4 dx, dz, rsnap, dsnap
CHARACTER*40 text
OPEN(10,form='formatted',file='snapshot.H')
text=' in="snapshot.H@"'
WRITE(10,*) text
text=' title="ELA2D snapshot data"'
WRITE(10,*) text
IF (nz.lt.10) THEN
WRITE(10,'(A4,I1)') ' n1=',nz
ELSEIF (nz.lt.100) THEN
WRITE(10,'(A4,I2)') ' n1=',nz
ELSEIF (nz.lt.1000) THEN
WRITE(10,'(A4,I3)') ' n1=',nz
ELSEIF (nz.lt.10000) THEN
WRITE(10,'(A4,I4)') ' n1=',nz
ENDIF
IF (nx.lt.10) THEN
WRITE(10,'(A4,I1)') ' n2=',nx
ELSEIF (nx.lt.100) THEN
WRITE(10,'(A4,I2)') ' n2=',nx
ELSEIF (nx.lt.1000) THEN
WRITE(10,'(A4,I3)') ' n2=',nx
ELSEIF (nx.lt.10000) THEN
WRITE(10,'(A4,I4)') ' n2=',nx
ENDIF
IF (nsnap.eq.0) THEN
WRITE(10,'(A4,I1)') ' n3=',0
ELSEIF (nsnap.lt.10) THEN
WRITE(10,'(A4,I1)') ' n3=',nsnap
ELSEIF (nsnap.lt.100) THEN
WRITE(10,'(A4,I2)') ' n3=',nsnap
ELSEIF (nsnap.lt.1000) THEN
WRITE(10,'(A4,I3)') ' n3=',nsnap
ENDIF
WRITE(10,*) ' esize=4'
WRITE(10,*) ' o1=0'
WRITE(10,*) ' o2=0'
IF (rsnap.lt.1.0) THEN
WRITE(10,'(A4,F6.5)') ' d3=',rsnap
ELSEIF (rsnap.lt.10.0) THEN
WRITE(10,'(A4,F6.4)') ' d3=',rsnap
ELSEIF (rsnap.lt.100.0) THEN
WRITE(10,'(A4,F6.3)') ' d3=',rsnap
ENDIF
IF (dsnap.lt.1.0) THEN
WRITE(10,'(A4,F8.7)') ' o3=',dsnap
ELSEIF (dsnap.lt.10.0) THEN
WRITE(10,'(A4,F8.6)') ' o3=',dsnap
ELSEIF (dsnap.lt.100.0) THEN
WRITE(10,'(A4,F8.5)') ' o3=',dsnap
ENDIF
IF (dx.lt.1.0) THEN
WRITE(10,'(A4,F8.7)') ' d2=',dx
ELSEIF (dx.lt.10.0) THEN
WRITE(10,'(A4,F8.6)') ' d2=',dx
ELSEIF (dx.lt.100.0) THEN
WRITE(10,'(A4,F8.5)') ' d2=',dx
ENDIF
IF (dz.lt.1.0) THEN
WRITE(10,'(A4,F6.5)') ' d1=',dz
ELSEIF (dz.lt.10.0) THEN
WRITE(10,'(A4,F6.4)') ' d1=',dz
ELSEIF (dz.lt.100.0) THEN
WRITE(10,'(A4,F6.3)') ' d1=',dz
ENDIF
WRITE(10,*) ' label1="m"'
WRITE(10,*) ' label2="m"'
WRITE(10,*) ' label3="sec"'
WRITE(10,*) ' data_format="xdr_float"'
CLOSE(10)
RETURN
END
English
