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

proceeding

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: 2274
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent a477349e
c this is <gresynoise.f>
c ----------------------------------------------------------------------------
c ($Id: gresynoise.f,v 1.1 2007-05-09 14:28:35 tforb Exp $)
c ($Id: gresynoise.f,v 1.2 2007-05-10 08:04:43 tforb Exp $)
c
c Copyright (c) 2007 by Thomas Forbriger (BFO Schiltach)
c
......@@ -10,6 +10,7 @@ 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 10/05/2007 V1.1 introduced randomization
c
c ============================================================================
c
......@@ -20,7 +21,7 @@ c
&'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 $')
&'$Id: gresynoise.f,v 1.2 2007-05-10 08:04:43 tforb Exp $')
c
c dimensions
integer maxtr, maxsamp, maxom, maxu
......@@ -36,6 +37,7 @@ c receiver specifications
real*8 rcvvred, rcvtli, rcvtre
character*80 rcvtext
real*8 phi(maxtr)
real scalefac(maxtr)
c seismogram parameters
integer pown
character*80 Zseisname,rcvname,Rseisname
......@@ -63,20 +65,21 @@ c seismogram trace to sff file
real dt
c commandline
integer maxopt, lastarg, iargc
parameter(maxopt=9)
parameter(maxopt=10)
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
logical suppress, optrndscale
real lambdalim, tapfrac, scalexp
c taper
real tf_costap
integer ltap, rtap
c source response
complex response(maxom)
real phase(maxom)
complex ime
parameter(ime=(0.,1.))
c any
integer lu, i, iargc, ierr, io, iu
real*8 pi2, arg, du, dom, scal, pi
......@@ -85,9 +88,10 @@ 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-/
data optid/2h-d,2h-l,2h-o,2h-t,2h-i,2h-S,2h-1,2h-2,2h-e,2h-r/
data opthasarg/.FALSE.,.TRUE.,.FALSE.,.TRUE.,4*.FALSE.,
& .TRUE.,.FALSE./
data optarg/1h-,2h1.,1h-,3h10.,4*1h-,2h0.,1h-/
c------------------------------------------------------------------------------
c basic information
c
......@@ -97,8 +101,8 @@ c
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 *,' [-t frac] [-e exp] [-r]'
print *,' [-S] [-1] [-2] [-i]'
print *,'or gresynoise -help'
if (greenname(1:5).ne.'-help')
& stop 'ERROR: wrong number of arguments'
......@@ -120,6 +124,11 @@ c
print *,' minimum wavelength in meters'
print *,' (default: ',optarg(2)(1:4),')'
print *,'-o overwrite output'
print *,'-i calculate impulse response'
print *,' do not randomize the signals''s phase'
print *,'-e exp scale seismograms with offset'
print *,' scaling factor: offset**exp'
print *,'-r additionally apply random scaling'
print *,'-t frac tapering fraction for slowness domain'
print *,' taper given in percent of full slowness'
print *,' range (default: ',optarg(4)(1:4),')'
......@@ -149,24 +158,20 @@ c
debug=optset(1)
optlambda=optset(2)
if (optlambda) then
read(optarg(2), *) lambdalim
read(optarg(2), *, err=99) lambdalim
endif
optnew=optset(3)
read(optarg(4), *)tapfrac
read(optarg(4), *, err=99) tapfrac
tapfrac=tapfrac/100.
optresponse=optset(5)
if (optresponse) stop 'ERROR: response file not supported'
respfile=optarg(5)
optresponse=.not.optset(5)
suppress=optset(6)
hankel1=optset(7)
hankel2=optset(8)
radial=optset(9)
if (radial) stop 'ERROR: radial option not supported'
read (optarg(9), *, err=99) scalexp
optscale=optset(10)
c
pi2=8.d0*datan(1.d0)
pi2=2.d0*pi
c----------------------------------------------------------------------
if (radial) print *,'calculate radial component'
c
c read configuration and greens function
c
call refmet_rrcv(rcvname, rcvtext,
......@@ -245,7 +250,9 @@ c
c scale down epicentral distances to meters
do i=1,ntr
r(i)=r(i)*1000.
scalefac(i)=1.
enddo
if (optscale) call tf_gsl_rnd_gaussian(scalefac, ntr)
c check stepwidth
dom=Zom(2)-Zom(1)
do io=2,Znom
......@@ -357,9 +364,10 @@ c----------------------------------------------------------------------
enddo
c apply response
if (optresponse) then
call tf_gsl_rng_uniform(phase, Znom)
do io=1,Znom
Zsdata(io)=Zsdata(io)*response(io)
Rsdata(io)=Rsdata(io)*response(io)
Zsdata(io)=Zsdata(io)*exp(ime*pi2*phase(io))
Rsdata(io)=Rsdata(io)*exp(ime*pi2*phase(io))
enddo
endif
do io=1,Znom-1
......@@ -372,8 +380,9 @@ c transform to time domain
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'
scal=1.d0/(sqrt(float(nsamp))*dt)*(r(i)**scalexp)
scal=scal*scalefac(i)
if (debug) print *,'stack traces'
do io=1,nsamp
Zfdata(io)=Zfdata(io)+real(Zsdata(io))*scal
if ((debug).and.
......
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