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

corrections

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: 2290
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 8eb18570
c this is <gresynoise.f>
c ----------------------------------------------------------------------------
c ($Id: gresynoise.f,v 1.11 2007-05-24 14:40:32 tforb Exp $)
c ($Id: gresynoise.f,v 1.12 2007-05-24 15:12:41 tforb Exp $)
c
c Copyright (c) 2007 by Thomas Forbriger (BFO Schiltach)
c
......@@ -15,6 +15,8 @@ c 11/05/2007 V1.2 allow selection of number of output samples
c 15/05/2007 V1.3 reordered loops; phase is of correct type now
c 21/05/2007 V1.4 provide longer seismograms by keep multiple copies of
c noisymized Fourier coefficients
c 24/05/2007 V1.5 write north and east component
c corrected buffer indexing
c
c ============================================================================
c
......@@ -45,10 +47,10 @@ c
c
character*(*) version
parameter(version=
&'GRESYNOISE V1.4 calculate noise seismograms')
&'GRESYNOISE V1.5 calculate noise seismograms')
character*(*) GRESYNOISE_CVS_ID
parameter(GRESYNOISE_CVS_ID=
&'$Id: gresynoise.f,v 1.11 2007-05-24 14:40:32 tforb Exp $')
&'$Id: gresynoise.f,v 1.12 2007-05-24 15:12:41 tforb Exp $')
c
c dimensions
integer maxtr, maxsetsamp, maxom, maxu, maxset, maxsamp
......@@ -93,6 +95,7 @@ c srce line
c info line
integer nstack
c seismogram trace to sff file
character cstation*5, cauxid*4, cinstype*6, cchannel*3
character*132 wid2line
logical last
real dt
......@@ -351,15 +354,22 @@ c
& 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
limnsetsamp=limnsetsamp/2
limnsetsamp=limnsetsamp
limnsetsamp=min(maxsetsamp/2,limnsetsamp)
limnsetsamp=max(Znom,limnsetsamp)
limnsetsamp=max(2*Znom,limnsetsamp)
limnsetsamp=limnsetsamp
pown=0
1 continue
pown=pown+1
nsetsamp=2**pown
if (nsetsamp.lt.limnsetsamp) goto 1
if (nsetsamp.lt.(2*Znom)) stop 'ERROR: cannot use enough samples'
if (nsetsamp.lt.(2*Znom)) then
print *,' limnsetsamp:',limnsetsamp
print *,' nsetsamp: ',nsetsamp
print *,' maxsetsamp: ',maxsetsamp
print *,' Znom: ',Znom
stop 'ERROR: cannot use enough samples'
endif
c calculate sampling interval
dt=pi2/(nsetsamp*dom)
nsamp=(nset+1)*nsetsamp/2
......@@ -387,7 +397,7 @@ c----------------------------------------------------------------------
c frequency loop
do io=1,Znom
c set slowness taper
if (debug) print *,'enter frequency loop'
cc if (debug) print *,'enter frequency loop'
if (optlambda) then
if (Zom(io).lt.dom) then
rtap=Znu
......@@ -412,7 +422,7 @@ c set slowness taper
c
c start slowness loop (conditional: select base function)
c -------------------------------------------------------
if (debug) print *,'enter slowness loop'
cc if (debug) print *,'enter slowness loop'
if (hankel1) then
do iu=1,rtap
arg=max(Zslo(iu)*Zom(io)*r(i),1.d-100)
......@@ -524,7 +534,7 @@ c extract time series
if (is.eq.nset) rtap=nsetsamp
do io=1,nsetsamp
scal=tf_costap(io,1,ltap,rtap,nsetsamp)
i=io+is*(nsetsamp/2)
i=io+(is-1)*(nsetsamp/2)
Zfdata(i)=Zfdata(i)+real(Ztrans(io))*scal
if ((debug).and.
& (abs(imag(Ztrans(io))).gt.(0.01*abs(real(Ztrans(io))))))
......@@ -559,7 +569,10 @@ c
c3=0.
cs='C'
nstack=0
last=.true.
last=.false.
cstation='NSP'
cauxid='NSP'
cinstype='NSP'
c vertical component
c ------------------
if (optnew) call sff_New(lu, Zseisname, ierr)
......@@ -570,16 +583,41 @@ c ------------------
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)
cchannel='Z'
call sff_prepwid2(nsamp, 1./dt, cstation, 1999, 4, 18, 0, 0,
& cchannel, cauxid, cinstype,
& 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 ----------------
c north component
c ---------------
c writing the same trace twice with equivalence buffers requires copying
do i=1,nsamp
Zfdata(i)=Rfdata(i)
enddo
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
cchannel='N'
call sff_prepwid2(nsamp, 1./dt, cstation, 1999, 4, 18, 0, 0,
& cchannel, cauxid, cinstype,
& 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 *,'written'
if (ierr.ne.0) stop 'ERROR: writing trace'
c east component
c --------------
last=.true.
if (optnew) call sff_New(lu, Rseisname, ierr)
if (ierr.ne.0) stop 'ERROR: deleting seismogram file'
call sff_WOpenFS(lu, Rseisname,
......@@ -588,8 +626,10 @@ c ----------------
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)
cchannel='E'
call sff_prepwid2(nsamp, 1./dt, cstation, 1999, 4, 18, 0, 0,
& cchannel, cauxid, cinstype,
& 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)
......
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