Commit a477349e authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

prototype version compiles

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 2273
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent d163f128
......@@ -73,6 +73,10 @@ syg: syg.o
gcc -o $(LOCBINDIR)/$@ $< -lgrrefsub -lrefread \
-ltf $(F2CLIB) -L$(LOCLIBDIR)
gresynoise: gresynoise.o
gcc -o $(LOCBINDIR)/$@ $< -lgrrefsub -lrefread \
-ltf -lgsl -lgslcblas $(LIBSTUFF) $(F2CLIB) -L$(LOCLIBDIR)
gresy: gresy.o
gcc -o $(LOCBINDIR)/$@ $< -lgrrefsub -lrefread \
-ltf $(LIBSTUFF) $(F2CLIB) -L$(LOCLIBDIR)
......
c this is <gresynoise.f>
c ----------------------------------------------------------------------------
c ($Id: gresynoise.f,v 1.1 2007-05-09 14:28:35 tforb Exp $)
c
c Copyright (c) 2007 by Thomas Forbriger (BFO Schiltach)
c
c calculate noise seismograms
c
c this file is obtained and modified from gresy.f (09/05/2007)
c
c REVISIONS and CHANGES
c 09/05/2007 V1.0 Thomas Forbriger
c
c ============================================================================
c
program gresynoise
c
character*(*) version
parameter(version=
&'GRESYNOISE V1.0 calculate noise seismograms')
character*(*) GRESYNOISE_CVS_ID
parameter(GRESYNOISE_CVS_ID=
&'$Id: gresynoise.f,v 1.1 2007-05-09 14:28:35 tforb Exp $')
c
c dimensions
integer maxtr, maxsamp, maxom, maxu
parameter(maxu=1000, maxtr=10000, maxom=4100, maxsamp=maxom*2)
integer ntr, nsamp, Znom, Znu, Rnom, Rnu
c greens function
character*78 greenname, Rgreenname
complex Zgreen(maxom, maxu)
complex Rgreen(maxom, maxu)
real Zom(maxom), Rom(maxom)
real Zslo(maxu), Rslo(maxu)
c receiver specifications
real*8 rcvvred, rcvtli, rcvtre
character*80 rcvtext
real*8 phi(maxtr)
c seismogram parameters
integer pown
character*80 Zseisname,rcvname,Rseisname
real*8 r(maxtr)
complex*16 Zsdata(maxsamp)
real Zfdata(maxsamp)
integer Zidata(maxsamp)
equivalence(Zfdata, Zidata)
complex*16 Rsdata(maxsamp)
real Rfdata(maxsamp)
integer Ridata(maxsamp)
equivalence(Rfdata, Ridata)
c free block
integer maxfree, nfree
parameter(maxfree=5)
character*80 free(maxfree)
c srce line
character srctype*40, date*20, time*20, cs
real c1, c2, c3
c info line
integer nstack
c seismogram trace to sff file
character*132 wid2line
logical last
real dt
c commandline
integer maxopt, lastarg, iargc
parameter(maxopt=9)
character*2 optid(maxopt)
character*80 optarg(maxopt)
logical optset(maxopt), opthasarg(maxopt)
c options
logical debug, optlambda, optnew, optresponse, hankel1, hankel2
logical suppress, radial
real lambdalim, tapfrac
character*80 respfile
c taper
real tf_costap
integer ltap, rtap
c source response
complex response(maxom)
c any
integer lu, i, iargc, ierr, io, iu
real*8 pi2, arg, du, dom, scal, pi
parameter(lu=20,pi=3.141592653589793115997)
c functions
real*8 tf_dj0, tf_dy0
real*8 tf_dj1, tf_dy1
c here are the keys to our commandline options
data optid/2h-d,2h-l,2h-o,2h-t,2h-r,2h-S,2h-1,2h-2,2h-R/
data opthasarg/.FALSE.,.TRUE.,.FALSE.,2*.TRUE.,4*.FALSE./
data optarg/1h-,2h1.,1h-,3h10.,4hjunk,4*1h-/
c------------------------------------------------------------------------------
c basic information
c
greenname=' '
if (iargc().eq.1) call getarg(1, greenname)
if ((greenname(1:5).eq.'-help').or.(iargc().lt.5)) then
print *,version
print *,'Usage: gresynoise Zgreen Rgreen Zseis Rseis rcvfile'
print *,' [-l lambda] [-o]'
print *,' [-t frac] [-r respfile]'
print *,' [-S] [-1] [-2] [-R]'
print *,'or gresynoise -help'
if (greenname(1:5).ne.'-help')
& stop 'ERROR: wrong number of arguments'
print *,' '
print *,GRESYNOISE_CVS_ID
print *,' '
print *,'Calculates synthetic noise seismograms'
print *,'from greens function.'
print *,' '
print *,'Zgreen file containing expension coeffcients'
print *,' produced by ''syg'' or ''greda'' '
print *,' for the vertical component'
print *,'Rgreen file containing expension coeffcients'
print *,' for the radial component'
print *,'Zseis output file for vertical component'
print *,'Rseis output file for radial component'
print *,'rcvfile receiver definition (refmet format)'
print *,'-l lambda limits the used slowness range by a'
print *,' minimum wavelength in meters'
print *,' (default: ',optarg(2)(1:4),')'
print *,'-o overwrite output'
print *,'-t frac tapering fraction for slowness domain'
print *,' taper given in percent of full slowness'
print *,' range (default: ',optarg(4)(1:4),')'
print *,'-S suppress zero frequency and zero slowness'
print *,'-1 use Hankel 1 instead of Bessel'
print *,'-2 use Hankel 2 instead of Bessel'
print *,' '
print *,'Input file units:'
print *,' '
print *,' Slowness unit: s/m'
print *,' Frequency unit: angular frequency 1/s'
print *,' amplitude unit: m**3/s if spectrum represents'
print *,' displacement waveform in m'
print *,' '
call refmet_rcvinf
stop
endif
if (iargc().lt.3) stop 'ERROR: need more arguments'
c
call getarg(2, Rgreenname)
call getarg(3, Zseisname)
call getarg(4, Rseisname)
call getarg(5, rcvname)
c
call tf_cmdline(6, lastarg,
& maxopt, optid, optarg, optset, opthasarg)
debug=optset(1)
optlambda=optset(2)
if (optlambda) then
read(optarg(2), *) lambdalim
endif
optnew=optset(3)
read(optarg(4), *)tapfrac
tapfrac=tapfrac/100.
optresponse=optset(5)
if (optresponse) stop 'ERROR: response file not supported'
respfile=optarg(5)
suppress=optset(6)
hankel1=optset(7)
hankel2=optset(8)
radial=optset(9)
if (radial) stop 'ERROR: radial option not supported'
c
pi2=8.d0*datan(1.d0)
c----------------------------------------------------------------------
if (radial) print *,'calculate radial component'
c
c read configuration and greens function
c
call refmet_rrcv(rcvname, rcvtext,
& rcvvred, rcvtli, rcvtre, ntr, maxtr, r, phi, 0.d0,
& 0, 1, debug)
c
print *,'receiver configuration was taken from:'
print *,' ',rcvtext
if (abs(rcvvred).gt.0.)
& print *,'WARNING: traveltime reduction will be ignored'
if (abs(rcvtli).gt.0.)
& print *,'WARNING: left time shift will be ignored'
if (abs(rcvtre).gt.0.)
& print *,'WARNING: right time shift will be ignored'
c
call greenread(greenname, debug,
& maxu, maxom, Zslo, Zom,
& Znu, Znom, Zgreen)
c
print *,'some information on green file ',
& 'read for vertical component:'
print *,' name: ',greenname(1:index(greenname, ' ')-1)
print *,' ',Znu,' slownesses from ',
& Zslo(1),'s/m to ',Zslo(Znu),'s/m'
print *,' ',Znom,' frequencies from ',
& Zom(1)/pi2,'Hz to ',Zom(Znom)/pi2,'Hz'
print *,' green frequency interval: ',(Zom(2)-Zom(1))/pi2,' Hz'
c
call greenread(Rgreenname, debug,
& maxu, maxom, Rslo, Rom,
& Rnu, Rnom, Rgreen)
c
print *,'some information on green file ',
& 'read for radial component:'
print *,' name: ',Rgreenname(1:index(Rgreenname, ' ')-1)
print *,' ',Rnu,' slownesses from ',
& Rslo(1),'s/m to ',Rslo(Rnu),'s/m'
print *,' ',Rnom,' frequencies from ',
& Rom(1)/pi2,'Hz to ',Rom(Rnom)/pi2,'Hz'
print *,' green frequency interval: ',(Rom(2)-Rom(1))/pi2,' Hz'
c
c both sets of coefficients must be coherent, since we have to apply the same
c noise signal
if (Znu.ne.Rnu) stop 'ERROR: numbers of slowness values differ'
if (Znom.ne.Rnom) stop 'ERROR: numbers of frequencies differ'
do io=1,maxom
if (abs(Rom(io)-Zom(io)).gt.(1.e-5*Zom(i)))
& stop 'ERROR: inconsistent frequency'
enddo
do iu=1,maxu
if (abs(Rslo(iu)-Zslo(iu)).gt.(1.e-5*Zslo(iu)))
& stop 'ERROR: inconsistent slowness'
enddo
c
if (suppress) then
print *,'suppress zero frequency and zero slowness coefficients'
do io=1,maxom
Zgreen(io,1)=(0.,0.)
Rgreen(io,1)=(0.,0.)
Zsdata(io)=(0.,0.)
Zfdata(io)=0.
Rsdata(io)=(0.,0.)
Rfdata(io)=0.
enddo
do iu=1,maxu
Zgreen(1,iu)=(0.,0.)
Rgreen(1,iu)=(0.,0.)
enddo
endif
c
c
c----------------------------------------------------------------------
c
c build seismic traces
c
c scale down epicentral distances to meters
do i=1,ntr
r(i)=r(i)*1000.
enddo
c check stepwidth
dom=Zom(2)-Zom(1)
do io=2,Znom
if (abs(Zom(io)-Zom(io-1)-dom).gt.(0.01*dom))
& print *,'WARNING: dom varies by more than 1%'
enddo
c
du=Zslo(2)-Zslo(1)
do iu=2,Znu
if (abs(Zslo(iu)-Zslo(iu-1)-du).gt.(0.01*du))
& print *,'WARNING: du varies by more than 1%'
enddo
c
if (Zom(1).gt.(0.01*dom))
& stop 'ERROR: I assume the first frequency to be 0.Hz!'
c we have got all frequencies now and are going to fit them to
c a power of 2 number
pown=0
1 continue
pown=pown+1
nsamp=2**pown
if (nsamp.lt.maxom) goto 1
c calculate sampling interval
dt=pi2/(nsamp*dom)
print *,'traces will have ',nsamp,' samples at ',dt,
& 's sampling interval'
c
if (hankel1) then
print *,'I will use the first Hankel function'
elseif (hankel2) then
print *,'I will use the second Hankel function'
else
print *,'I will use the Bessel function'
endif
last=.false.
do i=1,ntr
print *,'working on trace ',i,' at ',r(i),'m'
if (i.eq.ntr) last=.true.
c calculate bessel integral
c----------------------------------------------------------------------
c frequency loop
do io=1,Znom
c set slowness taper
if (debug) print *,'enter frequency loop'
if (optlambda) then
if (Zom(io).lt.dom) then
rtap=Znu
ltap=Znu
else
rtap=min(int((pi2/Zom(io)/lambdalim/du+1.)*
& (1.+tapfrac*0.5)),Znu)
ltap=min(int((pi2/Zom(io)/lambdalim/du+1.)*
& (1.-tapfrac*0.5)),Znu)
endif
ltap=min(ltap,int(rtap*(1.-tapfrac)))
else
rtap=Znu
ltap=int(rtap*(1.-tapfrac))
endif
if (i.eq.1) then
print *,'u-taper for ',Zom(io)/pi2,'Hz: ',ltap,'-',rtap,
& ' ',Zslo(ltap),'s/m-',Zslo(rtap),'s/m'
endif
Zsdata(io)=(0.d0,0.d0)
Rsdata(io)=(0.d0,0.d0)
if (debug) print *,'enter slowness loop'
if (hankel1) then
do iu=1,rtap
arg=max(Zslo(iu)*Zom(io)*r(i),1.d-100)
Rsdata(io)=Rsdata(io)+Rgreen(io,iu)*
& (tf_dj1(arg)+(0.d0,1.d0)*tf_dy1(arg))*
& du*tf_costap(iu,0,0,ltap,rtap)*Rslo(iu)
Zsdata(io)=Zsdata(io)+Zgreen(io,iu)*
& (tf_dj0(arg)+(0.d0,1.d0)*tf_dy0(arg))*
& du*tf_costap(iu,0,0,ltap,rtap)*Zslo(iu)
enddo
Zsdata(io)=0.5d0*Zsdata(io)
Rsdata(io)=0.5d0*Rsdata(io)
elseif (hankel2) then
do iu=1,rtap
arg=max(Zslo(iu)*Zom(io)*r(i),1.d-100)
Rsdata(io)=Rsdata(io)+Rgreen(io,iu)*
& (tf_dj1(arg)-(0.d0,1.d0)*tf_dy1(arg))*
& du*tf_costap(iu,0,0,ltap,rtap)*Rslo(iu)
Zsdata(io)=Zsdata(io)+Zgreen(io,iu)*
& (tf_dj0(arg)-(0.d0,1.d0)*tf_dy0(arg))*
& du*tf_costap(iu,0,0,ltap,rtap)*Zslo(iu)
enddo
Zsdata(io)=0.5d0*Zsdata(io)
Rsdata(io)=0.5d0*Rsdata(io)
else
do iu=1,rtap
arg=Zslo(iu)*Zom(io)*r(i)
Rsdata(io)=Rsdata(io)+Rgreen(io,iu)*
& tf_dj1(arg)*
& du*tf_costap(iu,0,0,ltap,rtap)*Rslo(iu)
Zsdata(io)=Zsdata(io)+Zgreen(io,iu)*
& tf_dj0(arg)*
& du*tf_costap(iu,0,0,ltap,rtap)*Zslo(iu)
enddo
endif
enddo
c end of frequency loop for Z-component
c----------------------------------------------------------------------
if (debug) print *,'left loops'
do io=Znom+1,nsamp-Znom+1
Zsdata(io)=(0.d0,0.d0)
Rsdata(io)=(0.d0,0.d0)
enddo
c apply response
if (optresponse) then
do io=1,Znom
Zsdata(io)=Zsdata(io)*response(io)
Rsdata(io)=Rsdata(io)*response(io)
enddo
endif
do io=1,Znom-1
Zsdata(nsamp-io+1)=conjg(Zsdata(io+1))
Rsdata(nsamp-io+1)=conjg(Rsdata(io+1))
enddo
c transform to time domain
if (debug) print *,'go fork'
call tf_dfork(nsamp, Zsdata, 1.d0)
call tf_dfork(nsamp, Rsdata, 1.d0)
if (debug) print *,'done fork'
c move trace to writing array and scale to impulse response
scal=1.d0/(sqrt(float(nsamp))*dt)
if (debug) print *,'move array'
do io=1,nsamp
Zfdata(io)=Zfdata(io)+real(Zsdata(io))*scal
if ((debug).and.
& (abs(imag(Zsdata(io))).gt.(0.01*abs(real(Zsdata(io))))))
& print *,'WHERE does this imaginary part come from?'
Rfdata(io)=Rfdata(io)+real(Rsdata(io))*scal
if ((debug).and.
& (abs(imag(Rsdata(io))).gt.(0.01*abs(real(Rsdata(io))))))
& print *,'WHERE does this imaginary part come from?'
enddo
enddo
c----------------------------------------------------------------------
c
c open seismogram file
c
free(1)=version
free(2)='green file '//greenname(1:index(greenname, ' '))
free(3)='receiver configuration from '//
& rcvname(1:index(rcvname, ' '))
free(4)=rcvtext
nfree=4
srctype='implicit in green file'
cs='C'
c1=0.
c2=0.
c3=0.
date='990418'
time='000000.000'
c1=0.
c2=0.
c3=0.
cs='C'
nstack=0
last=.true.
c vertical component
c ------------------
if (optnew) call sff_New(lu, Zseisname, ierr)
if (ierr.ne.0) stop 'ERROR: deleting seismogram file'
call sff_WOpenFS(lu, Zseisname,
& free, nfree, srctype, cs, c1 ,c2 ,c3,
& date, time, ierr)
if (ierr.ne.0) stop 'ERROR: opening seismogram file'
c prepare wid2line
if (debug) print *,'go and write',nsamp
call sff_prepwid2(nsamp, 1./dt, 'NSP', 1999, 4, 18, 0, 0, 'NSP ',
& 'NSP ', 'NSP ', 0., -1., -1., -1., -1., wid2line, ierr)
if (ierr.ne.0) stop 'ERROR: preparing WID2 line'
call sff_WTraceI(lu, wid2line, nsamp, Zfdata, Zidata, last,
& cs, c1, c2, c3, nstack, ierr)
if (debug) print *,'did writing'
if (ierr.ne.0) stop 'ERROR: writing trace'
c radial component
c ----------------
if (optnew) call sff_New(lu, Rseisname, ierr)
if (ierr.ne.0) stop 'ERROR: deleting seismogram file'
call sff_WOpenFS(lu, Rseisname,
& free, nfree, srctype, cs, c1 ,c2 ,c3,
& date, time, ierr)
if (ierr.ne.0) stop 'ERROR: opening seismogram file'
c prepare wid2line
if (debug) print *,'go and write',nsamp
call sff_prepwid2(nsamp, 1./dt, 'NSP', 1999, 4, 18, 0, 0, 'NSP ',
& 'NSP ', 'NSP ', 0., -1., -1., -1., -1., wid2line, ierr)
if (ierr.ne.0) stop 'ERROR: preparing WID2 line'
call sff_WTraceI(lu, wid2line, nsamp, Rfdata, Ridata, last,
& cs, c1, c2, c3, nstack, ierr)
if (debug) print *,'did writing'
if (ierr.ne.0) stop 'ERROR: writing trace'
c
stop
99 stop 'ERROR: reading command line argument'
98 stop 'ERROR: reading command line argument - to few distances'
97 stop 'ERROR opening response file'
96 stop 'ERROR reading response file'
95 stop 'ERROR reading response file - unexpected end'
94 stop 'ERROR closing response file'
50 format(//////)
end
c
c======================================================================
c
c some subroutines
c
c----------------------------------------------------------------------
c
c this routine reads in the data
c
subroutine greenread(filename, debug,
& maxslo, maxfreq, slo, om,
& nslo, nom, green)
c
character filename*(*)
logical debug
integer maxslo, maxfreq
complex green(maxfreq, maxslo)
real om(maxfreq), slo(maxslo)
integer nslo, nom
integer inmagic, cpu, match
character*4 cmagic, incmagic
parameter(cmagic='1234')
equivalence(incmagic, inmagic)
integer lu
parameter(lu=20)
integer i,j
c
print *,'read green file ',filename(1:index(filename, ' '))
open(lu, file=filename, form='unformatted', status='old', err=99)
c check byte sex
read(lu, err=98, end=97) inmagic
call tf_bytesex(cmagic, inmagic, cpu, match)
if (cpu.eq.1) then
print *,'running on Intel...'
elseif (cpu.eq.2) then
print *,'running on Motorola...'
else
stop 'ERROR: unknown processor type...'
endif
if (match.eq.1) then
print *,'matching bytesex - good...'
elseif (match.eq.2) then
print *,'bytesex not matching - we will have to swap!'
stop 'ERROR: do not know how to do that...'
else
close(lu, err=96)
print *,'bytesex read is ',incmagic
stop 'ERROR: bytesex is unkown - oh oh...'
endif
read(lu, err=98, end=97) nom, nslo
if ((nom.gt.maxfreq).or.(nslo.gt.maxslo)) then
close(lu, err=96)
stop 'ERROR: data exceeds array bounds'
endif
c read frequencies and phase slowness values
c frequency is given as angular frequency in 1/s
c phase slowness is given in s/m
read(lu, err=98, end=97) (om(i), i=1,nom), (slo(i), i=1,nslo)
read(lu, err=98, end=97) ((green(i,j), i=1,nom), j=1,nslo)
close(lu, err=96)
print *,'green file read and closed'
return
99 stop 'ERROR: opening green file'
98 stop 'ERROR: reading green file'
97 stop 'ERROR: reading green file - unexpected end'
96 stop 'ERROR: closing green file'
end
c ----- END OF gresy.f -----
c
c ----- END OF gresynoise.f -----
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment