mirror of
https://github.com/HongjianFang/DSurfTomo.git
synced 2025-05-09 00:51:12 +08:00
Compare commits
22 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
a651faf927 | ||
|
|
ad99372e34 | ||
|
|
37eb3e8f4b | ||
|
|
8d6695f308 | ||
|
|
65c202802c | ||
|
|
809e164bcb | ||
|
|
a0fceeb5fb | ||
|
|
2fec5df95b | ||
|
|
3d66dbdd0f | ||
|
|
f35c3d77da | ||
|
|
ba0a5db4ab | ||
|
|
b25eaaa096 | ||
|
|
f88ac4d8cd | ||
|
|
e4d1c1549f | ||
|
|
d8424b27b5 | ||
|
|
74e80f4b74 | ||
|
|
b657c68d28 | ||
|
|
9aead30bf2 | ||
|
|
82ebf4283d | ||
|
|
404f55cd0f | ||
|
|
f2a0f8bc4d | ||
|
|
07cf4d199b |
7
.gitignore
vendored
Normal file
7
.gitignore
vendored
Normal file
@ -0,0 +1,7 @@
|
||||
*.o
|
||||
*.mod
|
||||
Applications
|
||||
example_*
|
||||
.gitignore
|
||||
temp
|
||||
DSurfTomo
|
||||
21
README.md
21
README.md
@ -1,15 +1,8 @@
|
||||
#DSurfTomo
|
||||
DSurfTomo is the surface wave inversion program which can directly invert surface wave
|
||||
dispersion data to 3D shear wave speed without the intermediate step of constructing
|
||||
the phase or group velocity maps.A fast march method (FMM) is used to compute, at each
|
||||
period, surface wave travel times and ray paths between sources and receivers. This
|
||||
avoids the assumption of great-circle propagation that is used in most surface wave
|
||||
tomographic studies, but which is not appropriate in complex media.
|
||||
For detail description of the method, please refer to:
|
||||
Fang, H., Yao, H., Zhang, H., Huang, Y. C., & van der Hilst, R. D. (2015).
|
||||
Direct inversion of surface wave dispersion for three-dimensional shallow
|
||||
crustal structure based on ray tracing: methodology and application. Geophysical Journal International, 201(3), 1251-1263.
|
||||
DSurfTomo is a surface wave tomography program which inverts surface wave dispersion data directly to 3D shear wavespeed models without the intermediate step of constructing the phase or group velocity maps.
|
||||
The fast marching method (FMM) (Rawlinson et al., 2004) is used to compute, at each period, surface wave travel times and ray paths between sources and receivers. This avoids the assumption of great-circle propagation that is used in most surface wave tomographic studies, but which is not appropriate in complex media.
|
||||
To show its usage, an example of Taipei Basin tomography is provided, including scripts that are used for data reformating and plotting results. Interested users are recommended to refer to Fang et al. (2015 , GJI) for the detail description of the method.
|
||||
|
||||
Fang, H., Yao, H., Zhang, H., Huang, Y. C., & van der Hilst, R. D. (2015). Direct inversion of surface wave dispersion for three-dimensional shallow crustal structure based on ray tracing: methodology and application. Geophysical Journal International, 201(3), 1251-1263.
|
||||
|
||||
Rawlinson, N. & Sambridge, M., 2004. Wave front evolution in strongly heterogeneous layered media using the fast marching method, Geophys. J. Int., 156(3), 631–647
|
||||
|
||||
This is the source code without users' manual and examples about how to use it. If you need a detail description and some examples,
|
||||
feel free to send an email to Hongjian Fang (fanghj@mail.ustc.edu.cn) with a brief introduction about yourself and why you want to
|
||||
use it. We will be happy to help.
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
11
configure
vendored
11
configure
vendored
@ -1,16 +1,17 @@
|
||||
#!/usr/bin/env python
|
||||
# this is a script to install surf_tomo and surf_tomo_syn in your system
|
||||
# written by Hongjian Fang(fanghj@mail.ustc.edu.cn)
|
||||
|
||||
import os
|
||||
if 'bin' in os.listdir('.'):
|
||||
print 'installation beginning'
|
||||
print ('Installation beginning')
|
||||
else:
|
||||
print 'installation beginning'
|
||||
print ('Installation beginning')
|
||||
os.mkdir('bin')
|
||||
os.chdir('src')
|
||||
os.system('make clean')
|
||||
os.system('make')
|
||||
os.system('cp DSurfTomo ../bin')
|
||||
print '--------------------------------------'
|
||||
print 'Finishing DSurfTomo compiling'
|
||||
print '--------------------------------------'
|
||||
print ('--------------------------------------')
|
||||
print ('DSurfTomo compiling Finished')
|
||||
print ('--------------------------------------')
|
||||
|
||||
BIN
doc/ManualDSurfTomoV1.4.pdf
Normal file
BIN
doc/ManualDSurfTomoV1.4.pdf
Normal file
Binary file not shown.
@ -6,7 +6,7 @@ surfdataTB.dat c: data file
|
||||
25.2 121.35 c: goxd gozd (upper left point,[lat,lon])
|
||||
0.015 0.017 c: dvxd dvzd (grid interval in lat and lon direction)
|
||||
20 c: max(sources, receivers)
|
||||
4.0 0.1 c: weight damp
|
||||
4.0 1.0 c: weight damp
|
||||
3 c: sablayers (for computing depth kernel, 2~5)
|
||||
0.5 2.8 c: minimum velocity, maximum velocity (a priori information)
|
||||
10 c: maximum iteration
|
||||
@ -18,4 +18,4 @@ surfdataTB.dat c: data file
|
||||
0 c: kmaxLg
|
||||
0 c: synthetic flag(0:real data,1:synthetic)
|
||||
0.02 c: noiselevel
|
||||
0.5 c: threshold
|
||||
3.0 c: threshold
|
||||
146
scripts/GenerateDSurfTomoInputFile.py
Normal file
146
scripts/GenerateDSurfTomoInputFile.py
Normal file
@ -0,0 +1,146 @@
|
||||
#/* -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
|
||||
#
|
||||
# File Name : GenerateDSurfTomoInputFile.py
|
||||
#
|
||||
# Purpose : Generate the Input File for DSurfTomo using the data themselves
|
||||
#
|
||||
# Creation Date : 05-07-2023
|
||||
#
|
||||
# Last Modified : Tue 18 Jul 2023 08:10:22 PM CST
|
||||
#
|
||||
# Created By : Hongjian Fang: fanghj1990@gmail.com
|
||||
#
|
||||
#_._._._._._._._._._._._._._._._._._._._._.*/
|
||||
|
||||
import os
|
||||
import numpy as np
|
||||
import glob
|
||||
import pandas as pd
|
||||
|
||||
#These are for generating the right format
|
||||
DataPath = 'TaipeiRawData'
|
||||
DataNames = 'CDisp*.dat'
|
||||
wavetype = 2
|
||||
phasetype = 0
|
||||
|
||||
# These following can be changed in DSurfTomo.in
|
||||
nz = 9
|
||||
sublayers = 3
|
||||
minvel = 0.5
|
||||
maxvel = 2.8
|
||||
noiselevel = 0.02
|
||||
datathreshold = 1.5
|
||||
DataFiles = glob.glob(DataPath+'/'+DataNames)
|
||||
|
||||
sta1latall = []
|
||||
sta1lonall = []
|
||||
sta2latall = []
|
||||
sta2lonall = []
|
||||
periodsall = []
|
||||
dispersionall = []
|
||||
stationpairID = []
|
||||
|
||||
for ifile in DataFiles:
|
||||
dispersion = np.loadtxt(ifile)
|
||||
pairID = ifile.split('/')[-1]
|
||||
sta1lon = dispersion[0,0]
|
||||
sta1lat = dispersion[0,1]
|
||||
sta2lon = dispersion[1,0]
|
||||
sta2lat = dispersion[1,1]
|
||||
nperiods,_ = dispersion.shape
|
||||
nperiods = nperiods - 2
|
||||
periods = np.zeros(nperiods,)
|
||||
disper = np.zeros(nperiods,)
|
||||
for ii in range(nperiods):
|
||||
periods[ii] = dispersion[ii+2,0]
|
||||
disper[ii] = dispersion[ii+2,1]
|
||||
|
||||
dispersionall = np.hstack([dispersionall,disper])
|
||||
periodsall = np.hstack([periodsall,periods])
|
||||
sta1latall = np.hstack([sta1latall,np.ones(nperiods,)*sta1lat])
|
||||
sta1lonall = np.hstack([sta1lonall,np.ones(nperiods,)*sta1lon])
|
||||
sta2latall = np.hstack([sta2latall,np.ones(nperiods,)*sta2lat])
|
||||
sta2lonall = np.hstack([sta2lonall,np.ones(nperiods,)*sta2lon])
|
||||
stationpairID = stationpairID + [pairID] * nperiods
|
||||
|
||||
dataall = pd.DataFrame({'sta1lat':sta1latall, 'sta1lon':sta1lonall, \
|
||||
'sta2lat':sta2latall, 'sta2lon':sta2lonall, \
|
||||
'periods':periodsall, 'dispersion': dispersionall, \
|
||||
'pairid':stationpairID})
|
||||
|
||||
dataall = dataall.sort_values(by = ['periods', 'pairid'])
|
||||
dataall = dataall.set_index('periods')
|
||||
with open('surfdata.dat','w') as fout:
|
||||
UniqPeriods = dataall.index.unique()
|
||||
for iperiod,period in enumerate(UniqPeriods):
|
||||
datasubperiod = dataall.loc[period]
|
||||
if isinstance(datasubperiod,pd.Series):
|
||||
continue
|
||||
datasubperiod = datasubperiod.reset_index().set_index('sta1lat')
|
||||
sta1lat = datasubperiod.index.unique()
|
||||
for ista1 in sta1lat:
|
||||
datasubstation = datasubperiod.loc[ista1]
|
||||
if isinstance(datasubstation,pd.DataFrame):
|
||||
fout.write(f'# {datasubstation.index[0]} {datasubstation["sta1lon"].iloc[0]} {iperiod+1} {wavetype} {phasetype}\n')
|
||||
for ista2 in range(len(datasubstation)):
|
||||
fout.write(f'{datasubstation["sta2lat"].iloc[ista2]} {datasubstation["sta2lon"].iloc[ista2]} {datasubstation["dispersion"].iloc[ista2]}\n')
|
||||
else:
|
||||
fout.write(f'# {datasubstation.name} {datasubstation["sta1lon"]} {iperiod+1} {wavetype} {phasetype}\n')
|
||||
fout.write(f'{datasubstation["sta2lat"]} {datasubstation["sta2lon"]} {datasubstation["dispersion"]}\n')
|
||||
|
||||
|
||||
print('Finish reformatting disperison data')
|
||||
|
||||
# Determin the grid interval, 1/3 of the smallest station distance
|
||||
|
||||
# Weight and Damp, you will have to run the code a couple of time to decide
|
||||
# The default values is 5.0 and 1.0, respectively
|
||||
weight = 1.0
|
||||
damp = 1.0
|
||||
|
||||
# Sparsity fraction, this parameter represent the sparsity of the sensitivity matrix and is set to be 0.2 for safe keeping, but could be as low as 0.02.it is not important as long as your code does not throw you memory erros
|
||||
sparsityfraction = 0.1
|
||||
|
||||
# Iteration number, the default is 10 time, most often you get stable results after 4-5 iterations.
|
||||
maxiteration = 10
|
||||
|
||||
|
||||
distall = np.sqrt((dataall['sta2lat']-dataall['sta1lat'])**2+(dataall['sta2lon']-dataall['sta1lon'])**2)
|
||||
mindist = distall.min()
|
||||
gridintval = int(np.ceil(mindist/3*1000))/1000
|
||||
# Determine the origin, 1 grids
|
||||
originLat = np.max([np.max(dataall.sta1lat.max()),np.max(dataall.sta2lat.max())]) + 1*gridintval
|
||||
originLon = np.min([np.min(dataall.sta1lon.min()),np.min(dataall.sta2lon.min())]) - 1*gridintval
|
||||
largestLat = np.max([dataall['sta1lat'].max()-dataall['sta1lat'].min(),dataall['sta2lat'].max()-dataall['sta2lat'].min()])
|
||||
largestLon = np.max([dataall['sta1lon'].max()-dataall['sta1lon'].min(),dataall['sta2lon'].max()-dataall['sta2lon'].min()])
|
||||
inputfile = open('DSurfTomo.in','w')
|
||||
dummy = 'cccc'
|
||||
inputfile.write(f'{dummy*10}\n')
|
||||
inputfile.write(f'{dummy*10}\n')
|
||||
inputfile.write(f'{dummy*10}\n')
|
||||
inputfile.write('surfdata.dat\n')
|
||||
nx = int(np.ceil(largestLat/gridintval))+5
|
||||
ny = int(np.ceil(largestLon/gridintval))+5
|
||||
nreceivers = len(dataall.sta1lat.unique())+1
|
||||
inputfile.write(f'{nx} {ny} {nz}\n')
|
||||
inputfile.write(f'{originLat:<7.3f} {originLon:7.3f}\n')
|
||||
inputfile.write(f'{gridintval:<7.3f} {gridintval:7.3f}\n')
|
||||
inputfile.write(f'{nreceivers}\n')
|
||||
inputfile.write(f'{weight} {damp}\n')
|
||||
inputfile.write(f'{sublayers}\n')
|
||||
inputfile.write(f'{minvel} {maxvel}\n')
|
||||
inputfile.write(f'{maxiteration}\n')
|
||||
inputfile.write(f'{sparsityfraction}\n')
|
||||
nperiods = len(UniqPeriods)
|
||||
inputfile.write(f'{nperiods}\n')
|
||||
#inputfile.write(f'{*UniqPeriods}\n')
|
||||
inputfile.write(' '.join([str(iperiod) for iperiod in UniqPeriods])+'\n')
|
||||
inputfile.write(f'0\n')
|
||||
inputfile.write(f'0\n')
|
||||
inputfile.write(f'0\n')
|
||||
inputfile.write(f'0\n')
|
||||
inputfile.write(f'{noiselevel}\n')
|
||||
inputfile.write(f'{datathreshold}\n')
|
||||
inputfile.close()
|
||||
|
||||
print('Finishing generating input file\n')
|
||||
@ -8,10 +8,11 @@ import numpy as np
|
||||
#start
|
||||
nx=18
|
||||
ny=18
|
||||
nz=8
|
||||
minvel=0.8
|
||||
velgrad=0.5
|
||||
dep1=np.array([0,0.2,0.4,0.6,0.8,1.1,1.4,1.8,2.5])
|
||||
#dep1=np.array([0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.1,1.3,1.5,1.8,2.1,2.5])
|
||||
nz=len(dep1)
|
||||
#end
|
||||
vs1=np.zeros(nz)
|
||||
mod=np.zeros((nz*ny,nx))
|
||||
@ -29,4 +30,4 @@ with open('MOD','w') as fp:
|
||||
fp.write('%7.3f' % mod[k*ny+j,i])
|
||||
fp.write('\n')
|
||||
for i in range(nz):
|
||||
print dep1[i],
|
||||
print (dep1[i]),
|
||||
|
||||
69
scripts/plotcross.gmt
Normal file
69
scripts/plotcross.gmt
Normal file
@ -0,0 +1,69 @@
|
||||
#!/bin/csh
|
||||
|
||||
#-----------------------------------------------------------
|
||||
# 2015-06-08 Hongjian Fang
|
||||
# step 1: get velicoty along the track from each depth-layer
|
||||
# step 2: combine tracks into 2D profile
|
||||
#-----------------------------------------------------------
|
||||
|
||||
# change the following parameters according to your case
|
||||
# start
|
||||
set inp3D = DSurfTomo.inMeasure.dat
|
||||
set pstart = 121.4/25.1 #start point for 2D profile
|
||||
set pend = 121.55/25.1 #end point for 2D profile
|
||||
set pstart2 = 121.5/25.0 #start point for 2D profile
|
||||
set pend2 = 121.5/25.14 #end point for 2D profile
|
||||
set ddint = 0.017/0.015 #grid interval
|
||||
# over
|
||||
set diInt = 0.5 #distance interval (km)
|
||||
set J = -JX6.0i/-1.5i #size for plot
|
||||
set cpt = slice.cpt
|
||||
set ps = figcross.ps
|
||||
gmt makecpt -Cjet -I -T0.5/2.5/0.1 > $cpt #velocity boundary
|
||||
|
||||
# start
|
||||
set pro2d = tmp.profile.xzv
|
||||
set RMAP = `gmt gmtinfo -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
||||
rm -fr $pro2d
|
||||
|
||||
foreach layer(`awk '{print $3}' $inp3D| uniq`)
|
||||
set llayer = `echo $layer | awk '{printf "%04d",$1*1000}'`
|
||||
awk '{if($3==j) print $1,$2,$4}' j=$layer $inp3D |gmt surface -Gtmp.grd $RMAP -I$ddint -T0
|
||||
# be careful -I: grid spacing for x/y
|
||||
gmt project -C$pstart -E$pend -G$diInt -Q > track.dat
|
||||
gmt grdtrack track.dat -Gtmp.grd | awk '{print $3,j,$4}' j=$layer > tmp.D$llayer
|
||||
cat tmp.D$llayer >> $pro2d
|
||||
end
|
||||
#
|
||||
set RMAP = `gmt gmtinfo -C $pro2d | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
||||
gmt surface $pro2d -Gtmp.grd $RMAP -I$ddint -T0
|
||||
#gmt grdfilter tmp.grd -D0 -Fg4 -Gtmpf.grd
|
||||
gmt grdimage tmp.grd $J $RMAP -Bf2.5a5:'Distance (km)':/f0.5a1:'Depth (km)':SenW -C$cpt -K -Y5.2i> $ps
|
||||
#gmt psscale -Cslice.cpt -Ba0.5f0.25:'Vs': -D2.5i/-0.7i/6.00/0.2h -O -K -P >> $ps
|
||||
|
||||
rm -fr tmp.D* tmp.grd tmp.profile.xzv track.dat
|
||||
|
||||
# start
|
||||
set pro2d = tmp.profile.xzv
|
||||
set RMAP = `gmt gmtinfo -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
||||
rm -fr $pro2d
|
||||
|
||||
foreach layer(`awk '{print $3}' $inp3D| uniq`)
|
||||
set llayer = `echo $layer | awk '{printf "%04d",$1*1000}'`
|
||||
awk '{if($3==j) print $1,$2,$4}' j=$layer $inp3D |gmt surface -Gtmp.grd $RMAP -I$ddint -T0
|
||||
# be careful -I: grid spacing for x/y
|
||||
gmt project -C$pstart2 -E$pend2 -G$diInt -Q > track.dat
|
||||
gmt grdtrack track.dat -Gtmp.grd | awk '{print $3,j,$4}' j=$layer > tmp.D$llayer
|
||||
cat tmp.D$llayer >> $pro2d
|
||||
end
|
||||
#
|
||||
set RMAP = `gmt gmtinfo -C $pro2d | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
||||
gmt surface $pro2d -Gtmp.grd $RMAP -I$ddint -T0
|
||||
#gmt grdfilter tmp.grd -D0 -Fg4 -Gtmpf.grd
|
||||
gmt grdimage tmp.grd $J $RMAP -Bf2.5a5:'Distance (km)':/f0.5a1:'Depth (km)':SenW -C$cpt -K -O -Y-2.0i>> $ps
|
||||
gmt psscale -Cslice.cpt -Ba0.5f0.25:'Vs': -D2.5i/-0.7i/6.00/0.2h -O -P >> $ps
|
||||
|
||||
rm -fr tmp.D* tmp.grd tmp.profile.xzv track.dat
|
||||
p
|
||||
ps2pdf $ps
|
||||
|
||||
@ -7,42 +7,42 @@
|
||||
# csh plotslice.gmt SurfTomo.in.tvMeasure.dat 0.2 0.4 0.8 1.4
|
||||
#-----------------------------------------------------------
|
||||
|
||||
gmtset ANOT_FONT_SIZE 6
|
||||
gmtset LABEL_FONT_SIZE 6
|
||||
gmtset FRAME_WIDTH=0.1c
|
||||
gmtset LABEL_OFFSET=0.1c
|
||||
gmtset LABEL_FONT_SIZE=5p
|
||||
gmtset TICK_LENGTH=0.1c
|
||||
gmtset TICK_PEN=0.3p
|
||||
#gmt set ANOT_FONT_SIZE 6
|
||||
#gmt set LABEL_FONT_SIZE 6
|
||||
#gmt set FRAME_WIDTH=0.1c
|
||||
#gmt set LABEL_OFFSET=0.1c
|
||||
#gmt set LABEL_FONT_SIZE=5p
|
||||
#gmt set TICK_LENGTH=0.1c
|
||||
#gmt set TICK_PEN=0.3p
|
||||
set inp3D = $1
|
||||
set ps = horizVsnew.ps
|
||||
set ps = horizVsnew.true.ps
|
||||
set J = -JM2i #size for plot
|
||||
set cpt = slice.cpt
|
||||
|
||||
# start
|
||||
makecpt -Cseis -T0.6/1.5/0.1 > $cpt #velocity boundary
|
||||
set RMAP = `minmax -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
||||
psbasemap $RMAP $J -Ba0.1f0.05WseN -P -K -Y5i >$ps
|
||||
awk '{if($3==depth1) print $1,$2,$4}' depth1=$2 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||
gmt makecpt -Cseis -D -T0.5/2.5/0.1 > $cpt #velocity boundary
|
||||
set RMAP = `gmt gmtinfo -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
|
||||
gmt psbasemap $RMAP $J -Ba0.1f0.05WseN -P -K -Y5i >$ps
|
||||
awk '{if($3==depth1) print $1,$2,$4}' depth1=$2 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||
gmt grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||
rm -rf tmp.grd
|
||||
psscale -Cslice.cpt -Ba0.1f0.05:'S velocity (km/s)': -D1.0i/-0.15i/3.00/0.2h -O -K -P >> $ps
|
||||
gmt psscale -Cslice.cpt -Ba0.1f0.05:'S velocity (km/s)': -D1.0i/-0.15i/3.00/0.2h -O -K -P >> $ps
|
||||
|
||||
makecpt -Cseis -T0.8/1.8/0.1 > $cpt #velocity boundary
|
||||
psbasemap $RMAP $J -Ba0.1f0.05NwsE -P -K -O -X2.3i >>$ps
|
||||
awk '{if($3==depth2) print $1,$2,$4}' depth2=$3 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||
psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.15i/3.00/0.2h -O -K -P >> $ps
|
||||
#gmt makecpt -Cseis -D -T0.8/1.8/0.1 > $cpt #velocity boundary
|
||||
gmt psbasemap $RMAP $J -Ba0.1f0.05NwsE -P -K -O -X2.3i >>$ps
|
||||
awk '{if($3==depth2) print $1,$2,$4}' depth2=$3 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||
gmt grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||
gmt psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.15i/3.00/0.2h -O -K -P >> $ps
|
||||
|
||||
makecpt -Cseis -T1.1/2.0/0.1 > $cpt #velocity boundary
|
||||
psbasemap $RMAP $J -Ba0.1f0.05WneS -P -K -O -Y-2.7i -X-2.3i >>$ps
|
||||
awk '{if($3==depth2) print $1,$2,$4}' depth2=$4 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||
psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h -O -K -P >> $ps
|
||||
#gmt makecpt -Cseis -D -T1.1/2.0/0.1 > $cpt #velocity boundary
|
||||
gmt psbasemap $RMAP $J -Ba0.1f0.05WneS -P -K -O -Y-2.7i -X-2.3i >>$ps
|
||||
awk '{if($3==depth2) print $1,$2,$4}' depth2=$4 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||
gmt grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||
gmt psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h -O -K -P >> $ps
|
||||
|
||||
makecpt -Cseis -T1.3/2.4/0.1 > $cpt #velocity boundary
|
||||
psbasemap $RMAP $J -Ba0.1f0.05SwnE -P -K -O -X2.3i >>$ps
|
||||
awk '{if($3==depth2) print $1,$2,$4}' depth2=$5 $inp3D|xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||
grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||
psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h -O -P >> $ps
|
||||
#gmt makecpt -Cseis -D -T1.3/2.4/0.1 > $cpt #velocity boundary
|
||||
gmt psbasemap $RMAP $J -Ba0.1f0.05SwnE -P -K -O -X2.3i >>$ps
|
||||
awk '{if($3==depth2) print $1,$2,$4}' depth2=$5 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd
|
||||
gmt grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
|
||||
gmt psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h -O -P >> $ps
|
||||
|
||||
|
||||
@ -940,7 +940,7 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
||||
goxdf,gozdf,dvxdf,dvzdf,kmaxRc,kmaxRg,kmaxLc,kmaxLg, &
|
||||
tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk, &
|
||||
scxf,sczf,rcxf,rczf,nrc1,nsrcsurf1,kmax,nsrcsurf,nrcf, &
|
||||
nar,writepath)
|
||||
nar)
|
||||
USE globalp
|
||||
USE traveltime
|
||||
IMPLICIT NONE
|
||||
@ -1027,11 +1027,10 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
||||
integer ii,jj,kk,nn,istep
|
||||
integer level,maxlevel,maxleveld,HorizonType,VerticalType,PorS
|
||||
real,parameter::ftol=1e-4
|
||||
integer writepath
|
||||
integer ig, igroup
|
||||
|
||||
gdx=5
|
||||
gdz=5
|
||||
gdx=8
|
||||
gdz=8
|
||||
asgr=1
|
||||
sgdl=8
|
||||
sgs=8
|
||||
@ -1099,6 +1098,7 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
||||
if(kmaxRc.gt.0) then
|
||||
iwave=2
|
||||
igr=0
|
||||
print*,'Rayleigh wave phase velocity depth kernel'
|
||||
call depthkernel(nx,ny,nz,vels,pvRc,sen_vsRc,sen_vpRc, &
|
||||
sen_rhoRc,iwave,igr,kmaxRc,tRc,depz,minthk)
|
||||
endif
|
||||
@ -1108,8 +1108,9 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
||||
igr=0
|
||||
! print*,kmax
|
||||
call caldespersion(nx,ny,nz,vels,pvRc, &
|
||||
iwave,igr,kmax,tRg,depz,minthk)
|
||||
iwave,igr,kmaxRg,tRg,depz,minthk)
|
||||
igr=1
|
||||
print*,'Rayleigh wave group velocity depth kernel'
|
||||
call depthkernel(nx,ny,nz,vels,pvRg,sen_vsRg,sen_vpRg, &
|
||||
sen_rhoRg,iwave,igr,kmaxRg,tRg,depz,minthk)
|
||||
endif
|
||||
@ -1125,7 +1126,7 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
||||
iwave=1
|
||||
igr=0
|
||||
call caldespersion(nx,ny,nz,vels,pvLc, &
|
||||
iwave,igr,kmax,tLg,depz,minthk)
|
||||
iwave,igr,kmaxLg,tLg,depz,minthk)
|
||||
igr=1
|
||||
call depthkernel(nx,ny,nz,vels,pvLg,sen_vsLg,sen_vpLg, &
|
||||
sen_rhoLg,iwave,igr,kmaxLg,tLg,depz,minthk)
|
||||
@ -1282,6 +1283,7 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
||||
ALLOCATE(ttnr(idm2,idm1))
|
||||
ALLOCATE(nstsr(idm2,idm1))
|
||||
ENDIF
|
||||
!ttnr(1:nnzb,1:nnxb)=ttn(1:nnzb,1:nnxb)
|
||||
ttnr=ttn
|
||||
nstsr=nsts
|
||||
ogx=vnl
|
||||
@ -1376,13 +1378,11 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
||||
! if required.
|
||||
!
|
||||
if (igrt(srcnum,knumi) == 0 .or. (ig == 2 .and. igrt(srcnum,knumi) == 1)) then
|
||||
! a little stupid, remember to change latter
|
||||
if (igrt(srcnum,knumi) == 1) then
|
||||
call gridder(velf0)
|
||||
endif
|
||||
count11=count11+1
|
||||
CALL rpaths(x,z,fdm,rcxf(istep,srcnum,knumi),rczf(istep,srcnum,knumi),writepath)
|
||||
CALL rpaths(x,z,fdm,rcxf(istep,srcnum,knumi),rczf(istep,srcnum,knumi))
|
||||
row(1:nparpi)=0.0
|
||||
! change the Brocher relationship if depth is larger than 35 km
|
||||
if (depz(nz-1)<35.0) then
|
||||
do jj=1,nvz
|
||||
do kk=1,nvx
|
||||
if(abs(fdm(jj,kk)).ge.ftol) then
|
||||
@ -1400,6 +1400,28 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
! use a different formulation between Vp, Vs, and Rho, according to Luosong's
|
||||
else
|
||||
do jj=1,nvz
|
||||
do kk=1,nvx
|
||||
if(abs(fdm(jj,kk)).ge.ftol) then
|
||||
coe_a=(2.2110-0.8984*2*vels(kk+1,jj+1,1:nz-1)+&
|
||||
0.2786*3*vels(kk+1,jj+1,1:nz-1)**2-0.02412*4*vels(kk+1,jj+1,1:nz-1)**3)
|
||||
vpft=0.9098 + 2.2110*vels(kk+1,jj+1,1:nz-1) - 0.8984*vels(kk+1,jj+1,1:nz-1)**2+ &
|
||||
0.2786*vels(kk+1,jj+1,1:nz-1)**3 - 0.02412*vels(kk+1,jj+1,1:nz-1)**4
|
||||
coe_rho=coe_a*(1.6612-0.4721*2*vpft+&
|
||||
0.0671*3*vpft**2-0.0043*4*vpft**3+&
|
||||
0.000106*5*vpft**4)
|
||||
|
||||
row((jj-1)*nvx+kk:(nz-2)*nvz*nvx+(jj-1)*nvx+kk:nvx*nvz)=&
|
||||
(sen_vp(jj*(nvx+2)+kk+1,knumi,1:nz-1)*coe_a+&
|
||||
sen_rho(jj*(nvx+2)+kk+1,knumi,1:nz-1)*coe_rho+&
|
||||
sen_vs(jj*(nvx+2)+kk+1,knumi,1:nz-1))*fdm(jj,kk)
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
do nn=1,nparpi
|
||||
if(abs(row(nn)).gt.ftol) then
|
||||
nar=nar+1
|
||||
@ -1746,7 +1768,7 @@ END SUBROUTINE srtimes
|
||||
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
!SUBROUTINE rpaths(wrgf,csid,cfd,scx,scz)
|
||||
!SUBROUTINE rpaths()
|
||||
SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz,writepath)
|
||||
SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz)
|
||||
USE globalp
|
||||
IMPLICIT NONE
|
||||
INTEGER, PARAMETER :: i5=SELECTED_REAL_KIND(5,10)
|
||||
@ -1767,7 +1789,6 @@ SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz,writepath)
|
||||
!fang!------------------------------------------------
|
||||
real fdm(0:nvz+1,0:nvx+1)
|
||||
REAL(KIND=i10) surfrcx,surfrcz
|
||||
integer writepath
|
||||
!fang!------------------------------------------------
|
||||
!
|
||||
! ipx,ipz = Coordinates of cell containing current point
|
||||
@ -2252,14 +2273,14 @@ SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz,writepath)
|
||||
! Write ray paths to output file
|
||||
!
|
||||
!fang! IF(wrgf.EQ.csid.OR.wrgf.LT.0)THEN
|
||||
if(writepath == 1) then
|
||||
WRITE(40,*)'#',nrp
|
||||
DO j=1,nrp
|
||||
rayx=(pi/2-rgx(j))*180.0/pi
|
||||
rayz=rgz(j)*180.0/pi
|
||||
WRITE(40,*)rayx,rayz
|
||||
ENDDO
|
||||
endif
|
||||
!if(writepath == 1) then
|
||||
! WRITE(40,*)'#',nrp
|
||||
! DO j=1,nrp
|
||||
! rayx=(pi/2-rgx(j))*180.0/pi
|
||||
! rayz=rgz(j)*180.0/pi
|
||||
! WRITE(40,*)rayx,rayz
|
||||
! ENDDO
|
||||
!endif
|
||||
!fang! ENDIF
|
||||
!
|
||||
! Write partial derivatives to output file
|
||||
|
||||
12
src/Makefile
12
src/Makefile
@ -1,17 +1,19 @@
|
||||
CMD = DSurfTomo
|
||||
FC = gfortran
|
||||
FFLAGS = -O3 -ffixed-line-length-none -ffloat-store\
|
||||
-W -fbounds-check -m64 -mcmodel=medium
|
||||
FFLAGS = -O -ffixed-line-length-none -ffloat-store\
|
||||
-fbounds-check -m64 -mcmodel=medium
|
||||
F90SRCS = lsmrDataModule.f90 lsmrblasInterface.f90\
|
||||
lsmrblas.f90 lsmrModule.f90 delsph.f90\
|
||||
aprod.f90 gaussian.f90 main.f90
|
||||
FSRCS = surfdisp96.f
|
||||
OBJS = $(F90SRCS:%.f90=%.o) $(FSRCS:%.f=%.o) CalSurfG.o
|
||||
aprod.f90 gaussian.f90 getpercentile.f90
|
||||
FSRCS = surfdisp96.f slarnv.f slaruv.f
|
||||
OBJS = $(F90SRCS:%.f90=%.o) $(FSRCS:%.f=%.o) CalSurfG.o main.o
|
||||
all:$(CMD)
|
||||
$(CMD):$(OBJS)
|
||||
$(FC) -fopenmp $^ -o $@
|
||||
CalSurfG.o:CalSurfG.f90
|
||||
$(FC) -fopenmp $(FFLAGS) -c $< -o $@
|
||||
main.o:main.f90
|
||||
$(FC) -fopenmp $(FFLAGS) -c $< -o $@
|
||||
%.o: %.f90
|
||||
$(FC) $(FFLAGS) -c $(@F:.o=.f90) -o $@
|
||||
%.o: %.f
|
||||
|
||||
51
src/getpercentile.f90
Normal file
51
src/getpercentile.f90
Normal file
@ -0,0 +1,51 @@
|
||||
SUBROUTINE getpercentile(N,array,q25,q75)
|
||||
real q25
|
||||
real q75
|
||||
integer N
|
||||
real,intent(in) :: array(*)
|
||||
integer idx
|
||||
|
||||
real RA(N)
|
||||
RA(1:N) = array(1:N)
|
||||
L=N/2+1
|
||||
IR=N
|
||||
!The index L will be decremented from its initial value during the
|
||||
!"hiring" (heap creation) phase. Once it reaches 1, the index IR
|
||||
!will be decremented from its initial value down to 1 during the
|
||||
!"retirement-and-promotion" (heap selection) phase.
|
||||
10 continue
|
||||
if(L > 1)then
|
||||
L=L-1
|
||||
RRA=RA(L)
|
||||
else
|
||||
RRA=RA(IR)
|
||||
RA(IR)=RA(1)
|
||||
IR=IR-1
|
||||
if(IR.eq.1)then
|
||||
RA(1)=RRA
|
||||
idx = int(0.25*N)
|
||||
q25 = RA(idx)
|
||||
idx = int(0.75*N)
|
||||
q75 = RA(idx)
|
||||
return
|
||||
end if
|
||||
end if
|
||||
I=L
|
||||
J=L+L
|
||||
20 if(J.le.IR)then
|
||||
if(J < IR)then
|
||||
if(RA(J) < RA(J+1)) J=J+1
|
||||
end if
|
||||
if(RRA < RA(J))then
|
||||
RA(I)=RA(J)
|
||||
I=J; J=J+J
|
||||
else
|
||||
J=IR+1
|
||||
end if
|
||||
goto 20
|
||||
end if
|
||||
RA(I)=RRA
|
||||
goto 10
|
||||
|
||||
END
|
||||
|
||||
308
src/main.f90
308
src/main.f90
@ -20,6 +20,7 @@
|
||||
program SurfTomo
|
||||
use lsmrModule, only:lsmr
|
||||
use lsmrblasInterface, only : dnrm2
|
||||
use omp_lib
|
||||
implicit none
|
||||
|
||||
! VARIABLE DEFINE
|
||||
@ -62,8 +63,9 @@ program SurfTomo
|
||||
integer,dimension(:),allocatable::nsrc1
|
||||
integer,dimension(:,:),allocatable::periods
|
||||
real,dimension(:),allocatable::rw
|
||||
integer,dimension(:),allocatable::iw,col
|
||||
real,dimension(:),allocatable::dv,norm
|
||||
integer,dimension(:),allocatable::iw,col,nrow
|
||||
real,dimension(:),allocatable::dv,norm,dvsub,dvstd,dvall
|
||||
! real,dimension(:),allocatable::dvall
|
||||
real,dimension(:,:,:),allocatable::vsf
|
||||
real,dimension(:,:,:),allocatable::vsftrue
|
||||
character strf
|
||||
@ -96,34 +98,22 @@ program SurfTomo
|
||||
real spfra
|
||||
real noiselevel
|
||||
integer ifsyn
|
||||
integer writepath
|
||||
real averdws
|
||||
real maxnorm
|
||||
real threshold,threshold0
|
||||
real threshold0,q25,q75
|
||||
|
||||
! FOR MODEL VARIATION
|
||||
!------------------------------------------------
|
||||
integer idx
|
||||
integer counte
|
||||
real stdvs
|
||||
! integer numrand
|
||||
! real,allocatable,dimension(:,:)::modstat
|
||||
! real,allocatable,dimension(:)::modsig
|
||||
real gaussian
|
||||
external gaussian
|
||||
integer modest
|
||||
counte=0
|
||||
|
||||
! OPEN FILES FIRST TO OUTPUT THE PROCESS
|
||||
open(34,file='IterVel.out')
|
||||
nout=36
|
||||
open(nout,file='lsmr.txt')
|
||||
!nout=36
|
||||
!open(nout,file='lsmr.txt')
|
||||
|
||||
! OUTPUT PROGRAM INFOMATION
|
||||
write(*,*)
|
||||
write(*,*),' DSurfTomo (v1.3)'
|
||||
write(*,*),'PLEASE contact Hongjain Fang &
|
||||
(fanghj@mail.ustc.edu.cn) if you find any bug'
|
||||
write(*,*) ' DSurfTomo (v1.4)'
|
||||
!write(*,*) 'PLEASE contact Hongjain Fang &
|
||||
! (fanghj@mail.ustc.edu.cn) if you find any bug'
|
||||
write(*,*) 'For bug report, PLEASE contact Hongjain Fang &
|
||||
(fanghj1990@gmail.com)'
|
||||
write(*,*)
|
||||
|
||||
! READ INPUT FILE
|
||||
@ -145,7 +135,7 @@ program SurfTomo
|
||||
read(10,'(a30)')dummy
|
||||
read(10,'(a30)')dummy
|
||||
read(10,'(a30)')dummy
|
||||
read(10,*)datafile
|
||||
read(10,*) datafile
|
||||
read(10,*) nx,ny,nz
|
||||
read(10,*) goxd,gozd
|
||||
read(10,*) dvxd,dvzd
|
||||
@ -157,22 +147,22 @@ program SurfTomo
|
||||
read(10,*) spfra
|
||||
read(10,*) kmaxRc
|
||||
write(*,*) 'model origin:latitude,longitue'
|
||||
write(*,'(2f10.4)') goxd,gozd
|
||||
write(*,'(2f10.5)') goxd,gozd
|
||||
write(*,*) 'grid spacing:latitude,longitue'
|
||||
write(*,'(2f10.4)') dvxd,dvzd
|
||||
write(*,'(2f10.5)') dvxd,dvzd
|
||||
write(*,*) 'model dimension:nx,ny,nz'
|
||||
write(*,'(3i5)') nx,ny,nz
|
||||
write(logfile,'(a,a)')trim(inputfile),'.log'
|
||||
open(66,file=logfile)
|
||||
write(66,*)
|
||||
write(66,*),' S U R F T O M O'
|
||||
write(66,*),'PLEASE contact Hongjain Fang &
|
||||
write(66,*)' S U R F T O M O'
|
||||
write(66,*)'PLEASE contact Hongjain Fang &
|
||||
(fanghj@mail.ustc.edu.cn) if you find any bug'
|
||||
write(66,*)
|
||||
write(66,*) 'model origin:latitude,longitue'
|
||||
write(66,'(2f10.4)') goxd,gozd
|
||||
write(66,'(2f10.5)') goxd,gozd
|
||||
write(66,*) 'grid spacing:latitude,longitue'
|
||||
write(66,'(2f10.4)') dvxd,dvzd
|
||||
write(66,'(2f10.5)') dvxd,dvzd
|
||||
write(66,*) 'model dimension:nx,ny,nz'
|
||||
write(66,'(3i5)') nx,ny,nz
|
||||
if(kmaxRc.gt.0)then
|
||||
@ -181,9 +171,9 @@ program SurfTomo
|
||||
if (checkstat > 0) stop 'error allocating RP'
|
||||
read(10,*)(tRc(i),i=1,kmaxRc)
|
||||
write(*,*)'Rayleigh wave phase velocity used,periods:(s)'
|
||||
write(*,'(50f7.1)')(tRc(i),i=1,kmaxRc)
|
||||
write(*,'(50f7.2)')(tRc(i),i=1,kmaxRc)
|
||||
write(66,*)'Rayleigh wave phase velocity used,periods:(s)'
|
||||
write(66,'(50f7.1)')(tRc(i),i=1,kmaxRc)
|
||||
write(66,'(50f7.2)')(tRc(i),i=1,kmaxRc)
|
||||
endif
|
||||
read(10,*)kmaxRg
|
||||
if(kmaxRg.gt.0)then
|
||||
@ -191,9 +181,9 @@ program SurfTomo
|
||||
if (checkstat > 0) stop 'error allocating RP'
|
||||
read(10,*)(tRg(i),i=1,kmaxRg)
|
||||
write(*,*)'Rayleigh wave group velocity used,periods:(s)'
|
||||
write(*,'(50f7.1)')(tRg(i),i=1,kmaxRg)
|
||||
write(*,'(50f7.2)')(tRg(i),i=1,kmaxRg)
|
||||
write(66,*)'Rayleigh wave group velocity used,periods:(s)'
|
||||
write(66,'(50f7.1)')(tRg(i),i=1,kmaxRg)
|
||||
write(66,'(50f7.2)')(tRg(i),i=1,kmaxRg)
|
||||
endif
|
||||
read(10,*)kmaxLc
|
||||
if(kmaxLc.gt.0)then
|
||||
@ -201,9 +191,9 @@ program SurfTomo
|
||||
if (checkstat > 0) stop 'error allocating RP'
|
||||
read(10,*)(tLc(i),i=1,kmaxLc)
|
||||
write(*,*)'Love wave phase velocity used,periods:(s)'
|
||||
write(*,'(50f7.1)')(tLc(i),i=1,kmaxLc)
|
||||
write(*,'(50f7.2)')(tLc(i),i=1,kmaxLc)
|
||||
write(66,*)'Love wave phase velocity used,periods:(s)'
|
||||
write(66,'(50f7.1)')(tLc(i),i=1,kmaxLc)
|
||||
write(66,'(50f7.2)')(tLc(i),i=1,kmaxLc)
|
||||
endif
|
||||
read(10,*)kmaxLg
|
||||
if(kmaxLg.gt.0)then
|
||||
@ -211,15 +201,14 @@ program SurfTomo
|
||||
if (checkstat > 0) stop 'error allocating RP'
|
||||
read(10,*)(tLg(i),i=1,kmaxLg)
|
||||
write(*,*)'Love wave group velocity used,periods:(s)'
|
||||
write(*,'(50f7.1)')(tLg(i),i=1,kmaxLg)
|
||||
write(*,'(50f7.2)')(tLg(i),i=1,kmaxLg)
|
||||
write(66,*)'Love wave group velocity used,periods:(s)'
|
||||
write(66,'(50f7.1)')(tLg(i),i=1,kmaxLg)
|
||||
write(66,'(50f7.2)')(tLg(i),i=1,kmaxLg)
|
||||
endif
|
||||
read(10,*)ifsyn
|
||||
read(10,*)noiselevel
|
||||
read(10,*) threshold0
|
||||
! read(10,*) modest
|
||||
! read(10,*) numrand
|
||||
|
||||
close(10)
|
||||
nrc=nsrc
|
||||
kmax=kmaxRc+kmaxRg+kmaxLc+kmaxLg
|
||||
@ -296,15 +285,13 @@ program SurfTomo
|
||||
close(87)
|
||||
allocate(depz(nz), stat=checkstat)
|
||||
maxnar = spfra*dall*nx*ny*nz!sparsity fraction
|
||||
if (maxnar<0) print*, 'number overflow, decrease your sparsefrac'
|
||||
maxvp = (nx-2)*(ny-2)*(nz-1)
|
||||
allocate(dv(maxvp), stat=checkstat)
|
||||
allocate(dv(maxvp),dvsub(maxvp),dvstd(maxvp),dvall(maxvp), stat=checkstat)
|
||||
! allocate(dvall(maxvp*nrealizations),stats=checkstat)
|
||||
allocate(norm(maxvp), stat=checkstat)
|
||||
allocate(vsf(nx,ny,nz), stat=checkstat)
|
||||
allocate(vsftrue(nx,ny,nz), stat=checkstat)
|
||||
! FOR MODEL VARIATION
|
||||
!------------------------------------------------
|
||||
! allocate(modstat(numrand,maxvp))
|
||||
! allocate(modsig(maxvp))
|
||||
|
||||
allocate(rw(maxnar), stat=checkstat)
|
||||
if(checkstat > 0)then
|
||||
@ -314,7 +301,7 @@ program SurfTomo
|
||||
if(checkstat > 0)then
|
||||
write(6,*)'error with allocate: integer iw'
|
||||
endif
|
||||
allocate(col(maxnar), stat=checkstat)
|
||||
allocate(col(maxnar),nrow(dall), stat=checkstat)
|
||||
if(checkstat > 0)then
|
||||
write(6,*)'error with allocate: integer iw'
|
||||
endif
|
||||
@ -322,7 +309,7 @@ program SurfTomo
|
||||
stat=checkstat)
|
||||
|
||||
! MEASUREMENTS STATISTICS AND READ INITIAL MODEL
|
||||
write(*,'(a,i7)') 'Number of all measurements',dall
|
||||
write(*,'(a,i7)') ' Number of all measurements',dall
|
||||
|
||||
open(10,file='MOD',status='old')
|
||||
read(10,*) (depz(i),i=1,nz)
|
||||
@ -333,13 +320,11 @@ program SurfTomo
|
||||
enddo
|
||||
close(10)
|
||||
write(*,*) 'grid points in depth direction:(km)'
|
||||
write(*,'(50f7.1)') depz
|
||||
|
||||
|
||||
write(*,'(50f7.2)') depz
|
||||
|
||||
! CHECKERBOARD TEST
|
||||
if (ifsyn == 1) then
|
||||
write(*,*) 'Checkerboard Resolution Test Begin'
|
||||
write(*,*) 'Synthetic Test Begin'
|
||||
vsftrue = vsf
|
||||
|
||||
open(11,file='MOD.true',status='old')
|
||||
@ -360,60 +345,35 @@ program SurfTomo
|
||||
|
||||
|
||||
! ITERATE UNTILL CONVERGE
|
||||
writepath = 0
|
||||
do iter = 1,maxiter
|
||||
iw = 0
|
||||
rw = 0.0
|
||||
col = 0
|
||||
|
||||
! COMPUTE SENSITIVITY MATRIX
|
||||
if (iter == maxiter) then
|
||||
writepath = 1
|
||||
open(40,file='raypath.out')
|
||||
endif
|
||||
write(*,*) 'computing sensitivity matrix...'
|
||||
call CalSurfG(nx,ny,nz,maxvp,vsf,iw,rw,col,dsyn,&
|
||||
goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
|
||||
tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
|
||||
scxf,sczf,rcxf,rczf,nrc1,nsrc1,kmax,&
|
||||
nsrc,nrc,nar,writepath)
|
||||
nsrc,nrc,nar)
|
||||
|
||||
do i = 1,dall
|
||||
cbst(i) = obst(i) - dsyn(i)
|
||||
enddo
|
||||
|
||||
! write out rw, iw and cbst for testing
|
||||
!open(44,file='iw.dat')
|
||||
!write(44,*) nar
|
||||
!do i = 2,nar+1
|
||||
!write(44,*) iw(i)
|
||||
!enddo
|
||||
!do i = 1,nar
|
||||
!write(44,*) col(i)
|
||||
!enddo
|
||||
!close(44)
|
||||
|
||||
!open(44,file='rw.dat')
|
||||
!do i = 1,nar
|
||||
!write(44,*) rw(i)
|
||||
!enddo
|
||||
!close(44)
|
||||
|
||||
!open(44,file='residal.dat')
|
||||
!do i = 1,dall
|
||||
!write(44,*) cbst(i)
|
||||
!enddo
|
||||
!close(44)
|
||||
|
||||
!threshold = threshold0+(maxiter/2-iter)/3*0.5
|
||||
call getpercentile(dall,cbst,q25,q75)
|
||||
datweight = 1.0
|
||||
do i = 1,dall
|
||||
! datweight(i) = 1.0
|
||||
! if(abs(cbst(i)) > threshold) then
|
||||
! datweight(i) = exp(-(abs(cbst(i))-threshold))
|
||||
! endif
|
||||
datweight(i) = 0.01+1.0/(1+0.05*exp(cbst(i)**2*threshold0))
|
||||
cbst(i) = cbst(i)*datweight(i)
|
||||
if (cbst(i)<q25*threshold0 .or. cbst(i)>q75*threshold0) then
|
||||
datweight(i) = 0.0
|
||||
cbst(i) = 0
|
||||
endif
|
||||
enddo
|
||||
! do i = 1,dall
|
||||
! datweight(i) = 0.01+1.0/(1+0.05*exp(cbst(i)**2*threshold0))
|
||||
! cbst(i) = cbst(i)*datweight(i)
|
||||
! enddo
|
||||
|
||||
do i = 1,nar
|
||||
rw(i) = rw(i)*datweight(iw(1+i))
|
||||
@ -431,7 +391,6 @@ program SurfTomo
|
||||
enddo
|
||||
averdws=averdws/maxvp
|
||||
write(66,*)'Maximum and Average DWS values:',maxnorm,averdws
|
||||
write(66,*)'Threshold is:',threshold
|
||||
|
||||
! WRITE OUT RESIDUAL FOR THE FIRST AND LAST ITERATION
|
||||
if(iter.eq.1) then
|
||||
@ -452,8 +411,6 @@ program SurfTomo
|
||||
endif
|
||||
|
||||
|
||||
! ADDING REGULARIZATION TERM
|
||||
!weight=dnrm2(dall,cbst,1)**2/dall*weight0
|
||||
weight=weight0
|
||||
nar_tmp=nar
|
||||
nars=0
|
||||
@ -514,10 +471,12 @@ program SurfTomo
|
||||
leniw = 2*nar+1
|
||||
lenrw = nar
|
||||
dv = 0
|
||||
atol = 1e-3
|
||||
btol = 1e-3
|
||||
conlim = 1200
|
||||
itnlim = 1000
|
||||
atol = 1e-6
|
||||
btol = 1e-6
|
||||
!atol = 1e-3/((dvxd+dvzd)*111.19/2.0*0.1) !1e-2
|
||||
!btol = 1e-3/(dvxd*nx*111.19/3.0)!1e-3
|
||||
conlim = 100
|
||||
itnlim = 400
|
||||
istop = 0
|
||||
anorm = 0.0
|
||||
acond = 0.0
|
||||
@ -528,29 +487,32 @@ program SurfTomo
|
||||
call LSMR(m, n, leniw, lenrw,iw,rw,cbst, damp,&
|
||||
atol, btol, conlim, itnlim, localSize, nout,&
|
||||
dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
|
||||
if(istop==3) print*,'istop = 3, large condition number'
|
||||
|
||||
|
||||
do i =1,dall
|
||||
cbst(i)=cbst(i)/datweight(i)
|
||||
enddo
|
||||
|
||||
mean = sum(cbst(1:dall))/dall
|
||||
std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
|
||||
write(*,'(i2,a)'),iter,'th iteration...'
|
||||
write(*,'(a,f7.3)'),'weight is:',weight
|
||||
write(*,'(a,f8.1,a,f8.2,a,f8.3)'),'mean,std_devs and rms of &
|
||||
residual: ',mean*1000,'ms ',1000*std_devs,'ms ',&
|
||||
write(*,'(i2,a)')iter,'th iteration...'
|
||||
! write(*,'(a,f7.3)')'weight is:',weight
|
||||
write(*,'(a,f8.1,a,f8.2,a,f8.3)')' mean,std_devs and rms of &
|
||||
residual after weighting: ',mean*1000,'ms ',1000*std_devs,'ms ',&
|
||||
dnrm2(dall,cbst,1)/sqrt(real(dall))
|
||||
write(66,'(i2,a)'),iter,'th iteration...'
|
||||
write(66,'(a,f7.3)'),'weight is:',weight
|
||||
write(66,'(a,f8.1,a,f8.2,a,f8.3)'),'mean,std_devs and rms of &
|
||||
|
||||
!do i =1,dall
|
||||
!cbst(i)=cbst(i)/datweight(i)
|
||||
!enddo
|
||||
|
||||
!mean = sum(cbst(1:dall))/dall
|
||||
!std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
|
||||
!write(*,'(a,f8.1,a,f8.2,a,f8.3)')' residual before weighting: ',mean*1000,'ms ',1000*std_devs,'ms ',&
|
||||
! dnrm2(dall,cbst,1)/sqrt(real(dall))
|
||||
write(66,'(i2,a)')iter,'th iteration...'
|
||||
! write(66,'(a,f7.3)')'weight is:',weight
|
||||
write(66,'(a,f8.1,a,f8.2,a,f8.3)')'mean,std_devs and rms of &
|
||||
residual: ',mean*1000,'ms ',1000*std_devs,'ms ',&
|
||||
dnrm2(dall,cbst,1)/sqrt(real(dall))
|
||||
|
||||
write(*,'(a,2f7.4)'),'min and max velocity variation ',&
|
||||
write(*,'(a,2f7.4)')' min and max velocity variation ',&
|
||||
minval(dv),maxval(dv)
|
||||
write(66,'(a,2f7.4)'),'min and max velocity variation ',&
|
||||
write(66,'(a,2f7.4)')'min and max velocity variation ',&
|
||||
minval(dv),maxval(dv)
|
||||
|
||||
do k=1,nz-1
|
||||
@ -568,25 +530,13 @@ program SurfTomo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
write(34,*)'OUTPUT S VELOCITY AT ITERATION',iter
|
||||
do k=1,nz
|
||||
do j=1,ny
|
||||
write(34,'(100f7.3)') (vsf(i,j,k),i=1,nx)
|
||||
enddo
|
||||
enddo
|
||||
write(34,*)',OUTPUT DWS AT ITERATION',iter
|
||||
do k=1,nz-1
|
||||
do j=2,ny-1
|
||||
write(34,'(100f8.1)') (norm((k-1)*(ny-2)*(nx-2)+(j-2)*(nx-2)+i-1),i=2,nx-1)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
write(outmodel,'(a,a,i3.3)') trim(inputfile),'Measure.dat.iter',iter
|
||||
open(64,file=outmodel)
|
||||
do k=1,nz-1
|
||||
do j=1,ny-2
|
||||
do i=1,nx-2
|
||||
write(64,'(5f9.3)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i+1,j+1,k)
|
||||
write(64,'(5f10.5)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i+1,j+1,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -597,8 +547,8 @@ program SurfTomo
|
||||
|
||||
! OUTPUT THE VELOCITY MODEL
|
||||
|
||||
write(*,*),'Program finishes successfully'
|
||||
write(66,*),'Program finishes successfully'
|
||||
write(*,*)'Program finishes successfully'
|
||||
write(66,*)'Program finishes successfully'
|
||||
|
||||
if(ifsyn == 1) then
|
||||
open(65,file='Vs_model.real')
|
||||
@ -607,20 +557,20 @@ program SurfTomo
|
||||
do k=1,nz-1
|
||||
do j=1,ny-2
|
||||
do i=1,nx-2
|
||||
write(65,'(5f9.3)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsftrue(i+1,j+1,k)
|
||||
write(63,'(5f9.3)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i+1,j+1,k)
|
||||
write(65,'(5f10.5)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsftrue(i+1,j+1,k)
|
||||
write(63,'(5f10.5)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i+1,j+1,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(65)
|
||||
close(63)
|
||||
write(*,*),'Output True velocity model &
|
||||
write(*,*)'Output True velocity model &
|
||||
to Vs_model.real'
|
||||
write(*,*),'Output inverted shear velocity model &
|
||||
write(*,*)'Output inverted shear velocity model &
|
||||
to ',outsyn
|
||||
write(66,*),'Output True velocity model &
|
||||
write(66,*)'Output True velocity model &
|
||||
to Vs_model.real'
|
||||
write(66,*),'Output inverted shear velocity model &
|
||||
write(66,*)'Output inverted shear velocity model &
|
||||
to ',outsyn
|
||||
else
|
||||
write(outmodel,'(a,a)') trim(inputfile),'Measure.dat'
|
||||
@ -628,88 +578,22 @@ program SurfTomo
|
||||
do k=1,nz-1
|
||||
do j=1,ny-2
|
||||
do i=1,nx-2
|
||||
write(64,'(5f9.3)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i+1,j+1,k)
|
||||
write(64,'(5f10.5)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i+1,j+1,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(64)
|
||||
write(*,*),'Output inverted shear velocity model &
|
||||
write(*,*)'Output inverted shear velocity model &
|
||||
to ',outmodel
|
||||
write(66,*),'Output inverted shear velocity model &
|
||||
write(66,*)'Output inverted shear velocity model &
|
||||
to ',outmodel
|
||||
|
||||
|
||||
endif
|
||||
|
||||
close(34)
|
||||
close(40)
|
||||
close(nout) !close lsmr.txt
|
||||
!close(40)
|
||||
!close(nout) !close lsmr.txt
|
||||
close(66) !close surf_tomo.log
|
||||
|
||||
!! USE RANDOM MODEL TO OBTAIN THE MODEL VARIATION
|
||||
! !modest = 1
|
||||
! if (modest ==1) then
|
||||
!
|
||||
! write(*,*) 'model variation estimation begin...'
|
||||
! do iter = 1,numrand
|
||||
! call init_random_seed()
|
||||
! vsftrue=vsf
|
||||
! DO K=1,NZ-1
|
||||
! DO J=2,NY-1
|
||||
! DO I=2,NX-1
|
||||
! idx = (k-1)*(ny-2)*(nx-2)+(j-2)*(nx-2)+i-1
|
||||
! dv(idx) = 0.1/EXP(2*NORM(idx)/maxnorm)*gaussian()
|
||||
! VSFTRUE(I,J,K) = VSF(I,J,K)+dv(idx)
|
||||
! ENDDO
|
||||
! ENDDO
|
||||
! ENDDO
|
||||
! write(*,*),'maximum and minimum velocity variation',maxval(dv),minval(dv)
|
||||
!
|
||||
! call synthetic(nx,ny,nz,maxvp,vsftrue,dsyn,&
|
||||
! goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
|
||||
! tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
|
||||
! scxf,sczf,rcxf,rczf,nrc1,nsrc1,kmax,&
|
||||
! nsrc,nrc,0.0)
|
||||
!
|
||||
! do i = 1,dall
|
||||
! cbst(i) = obst(i) - dsyn(i)
|
||||
! enddo
|
||||
!
|
||||
! write(*,*), dnrm2(dall,cbst,1)/sqrt(real(dall)), 1.05*std_devs
|
||||
! if (dnrm2(dall,cbst,1)/sqrt(real(dall)) < 1.05*std_devs) then
|
||||
! counte = counte + 1
|
||||
! modstat(counte,:) = dv
|
||||
! endif
|
||||
!
|
||||
! enddo ! iteration for random models
|
||||
!
|
||||
! write(*,*),'number of of models satisfy requirements',counte
|
||||
! modsig = 1.0
|
||||
! if (counte>0) then
|
||||
! do i=1,maxvp
|
||||
! !statis
|
||||
! !mean = sum(cbst(1:dall))/dall
|
||||
! !std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
|
||||
! mean = sum(modstat(1:counte,i))/counte
|
||||
! stdvs = sqrt(sum(modstat(1:counte,i)**2)/counte-mean**2)
|
||||
! modsig(i) = stdvs
|
||||
! enddo
|
||||
! endif
|
||||
!
|
||||
! write(*,*),'write model variation to "model_variation.dat"'
|
||||
! open (64,file='model_variation.dat')
|
||||
! do k=1,nz-1
|
||||
! do j=1,ny-2
|
||||
! do i=1,nx-2
|
||||
! idx = (k-1)*(ny-2)*(nx-2)+(j-1)*(nx-2)+i
|
||||
! write(64,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),modsig(idx)
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! close(64)
|
||||
! write(*,*) 'finishing model variation estimation'
|
||||
! endif
|
||||
|
||||
|
||||
|
||||
deallocate(obst)
|
||||
deallocate(dsyn)
|
||||
deallocate(dist)
|
||||
@ -719,9 +603,9 @@ program SurfTomo
|
||||
deallocate(wavetype,igrt,nrc1)
|
||||
deallocate(nsrc1,periods)
|
||||
deallocate(rw)
|
||||
deallocate(iw,col)
|
||||
deallocate(iw,col,nrow)
|
||||
deallocate(cbst,wt,dtres,datweight)
|
||||
deallocate(dv)
|
||||
deallocate(dv,dvsub,dvstd,dvall)
|
||||
deallocate(norm)
|
||||
deallocate(vsf)
|
||||
deallocate(vsftrue)
|
||||
@ -740,19 +624,3 @@ program SurfTomo
|
||||
|
||||
end program
|
||||
|
||||
|
||||
|
||||
!-----------------------------------------------------------------------
|
||||
! Generate seed for random number generator of fortran
|
||||
! Note: only need to be called once inside one program
|
||||
!-----------------------------------------------------------------------
|
||||
subroutine init_random_seed()
|
||||
integer :: i,n,clock
|
||||
integer,dimension(:),allocatable :: seed
|
||||
call random_seed(size=n)
|
||||
allocate(seed(n))
|
||||
call system_clock(count=clock)
|
||||
seed=clock+37*(/(i-1,i=1,n)/)
|
||||
call random_seed(PUT=seed)
|
||||
deallocate(seed)
|
||||
end subroutine
|
||||
|
||||
178
src/slarnv.f
Normal file
178
src/slarnv.f
Normal file
@ -0,0 +1,178 @@
|
||||
*> \brief \b SLARNV returns a vector of random numbers from a uniform or normal distribution.
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
*
|
||||
* Online html documentation available at
|
||||
* http://www.netlib.org/lapack/explore-html/
|
||||
*
|
||||
*> \htmlonly
|
||||
*> Download SLARNV + dependencies
|
||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarnv.f">
|
||||
*> [TGZ]</a>
|
||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarnv.f">
|
||||
*> [ZIP]</a>
|
||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarnv.f">
|
||||
*> [TXT]</a>
|
||||
*> \endhtmlonly
|
||||
*
|
||||
* Definition:
|
||||
* ===========
|
||||
*
|
||||
* SUBROUTINE SLARNV( IDIST, ISEED, N, X )
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
* INTEGER IDIST, N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
* INTEGER ISEED( 4 )
|
||||
* REAL X( * )
|
||||
* ..
|
||||
*
|
||||
*
|
||||
*> \par Purpose:
|
||||
* =============
|
||||
*>
|
||||
*> \verbatim
|
||||
*>
|
||||
*> SLARNV returns a vector of n random real numbers from a uniform or
|
||||
*> normal distribution.
|
||||
*> \endverbatim
|
||||
*
|
||||
* Arguments:
|
||||
* ==========
|
||||
*
|
||||
*> \param[in] IDIST
|
||||
*> \verbatim
|
||||
*> IDIST is INTEGER
|
||||
*> Specifies the distribution of the random numbers:
|
||||
*> = 1: uniform (0,1)
|
||||
*> = 2: uniform (-1,1)
|
||||
*> = 3: normal (0,1)
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[in,out] ISEED
|
||||
*> \verbatim
|
||||
*> ISEED is INTEGER array, dimension (4)
|
||||
*> On entry, the seed of the random number generator; the array
|
||||
*> elements must be between 0 and 4095, and ISEED(4) must be
|
||||
*> odd.
|
||||
*> On exit, the seed is updated.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[in] N
|
||||
*> \verbatim
|
||||
*> N is INTEGER
|
||||
*> The number of random numbers to be generated.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[out] X
|
||||
*> \verbatim
|
||||
*> X is REAL array, dimension (N)
|
||||
*> The generated random numbers.
|
||||
*> \endverbatim
|
||||
*
|
||||
* Authors:
|
||||
* ========
|
||||
*
|
||||
*> \author Univ. of Tennessee
|
||||
*> \author Univ. of California Berkeley
|
||||
*> \author Univ. of Colorado Denver
|
||||
*> \author NAG Ltd.
|
||||
*
|
||||
*> \date December 2016
|
||||
*
|
||||
*> \ingroup OTHERauxiliary
|
||||
*
|
||||
*> \par Further Details:
|
||||
* =====================
|
||||
*>
|
||||
*> \verbatim
|
||||
*>
|
||||
*> This routine calls the auxiliary routine SLARUV to generate random
|
||||
*> real numbers from a uniform (0,1) distribution, in batches of up to
|
||||
*> 128 using vectorisable code. The Box-Muller method is used to
|
||||
*> transform numbers from a uniform to a normal distribution.
|
||||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE SLARNV( IDIST, ISEED, N, X )
|
||||
*
|
||||
* -- LAPACK auxiliary routine (version 3.7.0) --
|
||||
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
* December 2016
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER IDIST, N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
INTEGER ISEED( 4 )
|
||||
REAL X( * )
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
REAL ONE, TWO
|
||||
PARAMETER ( ONE = 1.0E+0, TWO = 2.0E+0 )
|
||||
INTEGER LV
|
||||
PARAMETER ( LV = 128 )
|
||||
REAL TWOPI
|
||||
PARAMETER ( TWOPI = 6.2831853071795864769252867663E+0 )
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
INTEGER I, IL, IL2, IV
|
||||
* ..
|
||||
* .. Local Arrays ..
|
||||
REAL U( LV )
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC COS, LOG, MIN, SQRT
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL SLARUV
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
DO 40 IV = 1, N, LV / 2
|
||||
IL = MIN( LV / 2, N-IV+1 )
|
||||
IF( IDIST.EQ.3 ) THEN
|
||||
IL2 = 2*IL
|
||||
ELSE
|
||||
IL2 = IL
|
||||
END IF
|
||||
*
|
||||
* Call SLARUV to generate IL2 numbers from a uniform (0,1)
|
||||
* distribution (IL2 <= LV)
|
||||
*
|
||||
CALL SLARUV( ISEED, IL2, U )
|
||||
*
|
||||
IF( IDIST.EQ.1 ) THEN
|
||||
*
|
||||
* Copy generated numbers
|
||||
*
|
||||
DO 10 I = 1, IL
|
||||
X( IV+I-1 ) = U( I )
|
||||
10 CONTINUE
|
||||
ELSE IF( IDIST.EQ.2 ) THEN
|
||||
*
|
||||
* Convert generated numbers to uniform (-1,1) distribution
|
||||
*
|
||||
DO 20 I = 1, IL
|
||||
X( IV+I-1 ) = TWO*U( I ) - ONE
|
||||
20 CONTINUE
|
||||
ELSE IF( IDIST.EQ.3 ) THEN
|
||||
*
|
||||
* Convert generated numbers to normal (0,1) distribution
|
||||
*
|
||||
DO 30 I = 1, IL
|
||||
X( IV+I-1 ) = SQRT( -TWO*LOG( U( 2*I-1 ) ) )*
|
||||
$ COS( TWOPI*U( 2*I ) )
|
||||
30 CONTINUE
|
||||
END IF
|
||||
40 CONTINUE
|
||||
RETURN
|
||||
*
|
||||
* End of SLARNV
|
||||
*
|
||||
END
|
||||
447
src/slaruv.f
Normal file
447
src/slaruv.f
Normal file
@ -0,0 +1,447 @@
|
||||
*> \brief \b SLARUV returns a vector of n random real numbers from a uniform distribution.
|
||||
*
|
||||
* =========== DOCUMENTATION ===========
|
||||
*
|
||||
* Online html documentation available at
|
||||
* http://www.netlib.org/lapack/explore-html/
|
||||
*
|
||||
*> \htmlonly
|
||||
*> Download SLARUV + dependencies
|
||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaruv.f">
|
||||
*> [TGZ]</a>
|
||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaruv.f">
|
||||
*> [ZIP]</a>
|
||||
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaruv.f">
|
||||
*> [TXT]</a>
|
||||
*> \endhtmlonly
|
||||
*
|
||||
* Definition:
|
||||
* ===========
|
||||
*
|
||||
* SUBROUTINE SLARUV( ISEED, N, X )
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
* INTEGER N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
* INTEGER ISEED( 4 )
|
||||
* REAL X( N )
|
||||
* ..
|
||||
*
|
||||
*
|
||||
*> \par Purpose:
|
||||
* =============
|
||||
*>
|
||||
*> \verbatim
|
||||
*>
|
||||
*> SLARUV returns a vector of n random real numbers from a uniform (0,1)
|
||||
*> distribution (n <= 128).
|
||||
*>
|
||||
*> This is an auxiliary routine called by SLARNV and CLARNV.
|
||||
*> \endverbatim
|
||||
*
|
||||
* Arguments:
|
||||
* ==========
|
||||
*
|
||||
*> \param[in,out] ISEED
|
||||
*> \verbatim
|
||||
*> ISEED is INTEGER array, dimension (4)
|
||||
*> On entry, the seed of the random number generator; the array
|
||||
*> elements must be between 0 and 4095, and ISEED(4) must be
|
||||
*> odd.
|
||||
*> On exit, the seed is updated.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[in] N
|
||||
*> \verbatim
|
||||
*> N is INTEGER
|
||||
*> The number of random numbers to be generated. N <= 128.
|
||||
*> \endverbatim
|
||||
*>
|
||||
*> \param[out] X
|
||||
*> \verbatim
|
||||
*> X is REAL array, dimension (N)
|
||||
*> The generated random numbers.
|
||||
*> \endverbatim
|
||||
*
|
||||
* Authors:
|
||||
* ========
|
||||
*
|
||||
*> \author Univ. of Tennessee
|
||||
*> \author Univ. of California Berkeley
|
||||
*> \author Univ. of Colorado Denver
|
||||
*> \author NAG Ltd.
|
||||
*
|
||||
*> \date December 2016
|
||||
*
|
||||
*> \ingroup OTHERauxiliary
|
||||
*
|
||||
*> \par Further Details:
|
||||
* =====================
|
||||
*>
|
||||
*> \verbatim
|
||||
*>
|
||||
*> This routine uses a multiplicative congruential method with modulus
|
||||
*> 2**48 and multiplier 33952834046453 (see G.S.Fishman,
|
||||
*> 'Multiplicative congruential random number generators with modulus
|
||||
*> 2**b: an exhaustive analysis for b = 32 and a partial analysis for
|
||||
*> b = 48', Math. Comp. 189, pp 331-344, 1990).
|
||||
*>
|
||||
*> 48-bit integers are stored in 4 integer array elements with 12 bits
|
||||
*> per element. Hence the routine is portable across machines with
|
||||
*> integers of 32 bits or more.
|
||||
*> \endverbatim
|
||||
*>
|
||||
* =====================================================================
|
||||
SUBROUTINE SLARUV( ISEED, N, X )
|
||||
*
|
||||
* -- LAPACK auxiliary routine (version 3.7.0) --
|
||||
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
||||
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||||
* December 2016
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
INTEGER ISEED( 4 )
|
||||
REAL X( N )
|
||||
* ..
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
REAL ONE
|
||||
PARAMETER ( ONE = 1.0E0 )
|
||||
INTEGER LV, IPW2
|
||||
REAL R
|
||||
PARAMETER ( LV = 128, IPW2 = 4096, R = ONE / IPW2 )
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
INTEGER I, I1, I2, I3, I4, IT1, IT2, IT3, IT4, J
|
||||
* ..
|
||||
* .. Local Arrays ..
|
||||
INTEGER MM( LV, 4 )
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MIN, MOD, REAL
|
||||
* ..
|
||||
* .. Data statements ..
|
||||
DATA ( MM( 1, J ), J = 1, 4 ) / 494, 322, 2508,
|
||||
$ 2549 /
|
||||
DATA ( MM( 2, J ), J = 1, 4 ) / 2637, 789, 3754,
|
||||
$ 1145 /
|
||||
DATA ( MM( 3, J ), J = 1, 4 ) / 255, 1440, 1766,
|
||||
$ 2253 /
|
||||
DATA ( MM( 4, J ), J = 1, 4 ) / 2008, 752, 3572,
|
||||
$ 305 /
|
||||
DATA ( MM( 5, J ), J = 1, 4 ) / 1253, 2859, 2893,
|
||||
$ 3301 /
|
||||
DATA ( MM( 6, J ), J = 1, 4 ) / 3344, 123, 307,
|
||||
$ 1065 /
|
||||
DATA ( MM( 7, J ), J = 1, 4 ) / 4084, 1848, 1297,
|
||||
$ 3133 /
|
||||
DATA ( MM( 8, J ), J = 1, 4 ) / 1739, 643, 3966,
|
||||
$ 2913 /
|
||||
DATA ( MM( 9, J ), J = 1, 4 ) / 3143, 2405, 758,
|
||||
$ 3285 /
|
||||
DATA ( MM( 10, J ), J = 1, 4 ) / 3468, 2638, 2598,
|
||||
$ 1241 /
|
||||
DATA ( MM( 11, J ), J = 1, 4 ) / 688, 2344, 3406,
|
||||
$ 1197 /
|
||||
DATA ( MM( 12, J ), J = 1, 4 ) / 1657, 46, 2922,
|
||||
$ 3729 /
|
||||
DATA ( MM( 13, J ), J = 1, 4 ) / 1238, 3814, 1038,
|
||||
$ 2501 /
|
||||
DATA ( MM( 14, J ), J = 1, 4 ) / 3166, 913, 2934,
|
||||
$ 1673 /
|
||||
DATA ( MM( 15, J ), J = 1, 4 ) / 1292, 3649, 2091,
|
||||
$ 541 /
|
||||
DATA ( MM( 16, J ), J = 1, 4 ) / 3422, 339, 2451,
|
||||
$ 2753 /
|
||||
DATA ( MM( 17, J ), J = 1, 4 ) / 1270, 3808, 1580,
|
||||
$ 949 /
|
||||
DATA ( MM( 18, J ), J = 1, 4 ) / 2016, 822, 1958,
|
||||
$ 2361 /
|
||||
DATA ( MM( 19, J ), J = 1, 4 ) / 154, 2832, 2055,
|
||||
$ 1165 /
|
||||
DATA ( MM( 20, J ), J = 1, 4 ) / 2862, 3078, 1507,
|
||||
$ 4081 /
|
||||
DATA ( MM( 21, J ), J = 1, 4 ) / 697, 3633, 1078,
|
||||
$ 2725 /
|
||||
DATA ( MM( 22, J ), J = 1, 4 ) / 1706, 2970, 3273,
|
||||
$ 3305 /
|
||||
DATA ( MM( 23, J ), J = 1, 4 ) / 491, 637, 17,
|
||||
$ 3069 /
|
||||
DATA ( MM( 24, J ), J = 1, 4 ) / 931, 2249, 854,
|
||||
$ 3617 /
|
||||
DATA ( MM( 25, J ), J = 1, 4 ) / 1444, 2081, 2916,
|
||||
$ 3733 /
|
||||
DATA ( MM( 26, J ), J = 1, 4 ) / 444, 4019, 3971,
|
||||
$ 409 /
|
||||
DATA ( MM( 27, J ), J = 1, 4 ) / 3577, 1478, 2889,
|
||||
$ 2157 /
|
||||
DATA ( MM( 28, J ), J = 1, 4 ) / 3944, 242, 3831,
|
||||
$ 1361 /
|
||||
DATA ( MM( 29, J ), J = 1, 4 ) / 2184, 481, 2621,
|
||||
$ 3973 /
|
||||
DATA ( MM( 30, J ), J = 1, 4 ) / 1661, 2075, 1541,
|
||||
$ 1865 /
|
||||
DATA ( MM( 31, J ), J = 1, 4 ) / 3482, 4058, 893,
|
||||
$ 2525 /
|
||||
DATA ( MM( 32, J ), J = 1, 4 ) / 657, 622, 736,
|
||||
$ 1409 /
|
||||
DATA ( MM( 33, J ), J = 1, 4 ) / 3023, 3376, 3992,
|
||||
$ 3445 /
|
||||
DATA ( MM( 34, J ), J = 1, 4 ) / 3618, 812, 787,
|
||||
$ 3577 /
|
||||
DATA ( MM( 35, J ), J = 1, 4 ) / 1267, 234, 2125,
|
||||
$ 77 /
|
||||
DATA ( MM( 36, J ), J = 1, 4 ) / 1828, 641, 2364,
|
||||
$ 3761 /
|
||||
DATA ( MM( 37, J ), J = 1, 4 ) / 164, 4005, 2460,
|
||||
$ 2149 /
|
||||
DATA ( MM( 38, J ), J = 1, 4 ) / 3798, 1122, 257,
|
||||
$ 1449 /
|
||||
DATA ( MM( 39, J ), J = 1, 4 ) / 3087, 3135, 1574,
|
||||
$ 3005 /
|
||||
DATA ( MM( 40, J ), J = 1, 4 ) / 2400, 2640, 3912,
|
||||
$ 225 /
|
||||
DATA ( MM( 41, J ), J = 1, 4 ) / 2870, 2302, 1216,
|
||||
$ 85 /
|
||||
DATA ( MM( 42, J ), J = 1, 4 ) / 3876, 40, 3248,
|
||||
$ 3673 /
|
||||
DATA ( MM( 43, J ), J = 1, 4 ) / 1905, 1832, 3401,
|
||||
$ 3117 /
|
||||
DATA ( MM( 44, J ), J = 1, 4 ) / 1593, 2247, 2124,
|
||||
$ 3089 /
|
||||
DATA ( MM( 45, J ), J = 1, 4 ) / 1797, 2034, 2762,
|
||||
$ 1349 /
|
||||
DATA ( MM( 46, J ), J = 1, 4 ) / 1234, 2637, 149,
|
||||
$ 2057 /
|
||||
DATA ( MM( 47, J ), J = 1, 4 ) / 3460, 1287, 2245,
|
||||
$ 413 /
|
||||
DATA ( MM( 48, J ), J = 1, 4 ) / 328, 1691, 166,
|
||||
$ 65 /
|
||||
DATA ( MM( 49, J ), J = 1, 4 ) / 2861, 496, 466,
|
||||
$ 1845 /
|
||||
DATA ( MM( 50, J ), J = 1, 4 ) / 1950, 1597, 4018,
|
||||
$ 697 /
|
||||
DATA ( MM( 51, J ), J = 1, 4 ) / 617, 2394, 1399,
|
||||
$ 3085 /
|
||||
DATA ( MM( 52, J ), J = 1, 4 ) / 2070, 2584, 190,
|
||||
$ 3441 /
|
||||
DATA ( MM( 53, J ), J = 1, 4 ) / 3331, 1843, 2879,
|
||||
$ 1573 /
|
||||
DATA ( MM( 54, J ), J = 1, 4 ) / 769, 336, 153,
|
||||
$ 3689 /
|
||||
DATA ( MM( 55, J ), J = 1, 4 ) / 1558, 1472, 2320,
|
||||
$ 2941 /
|
||||
DATA ( MM( 56, J ), J = 1, 4 ) / 2412, 2407, 18,
|
||||
$ 929 /
|
||||
DATA ( MM( 57, J ), J = 1, 4 ) / 2800, 433, 712,
|
||||
$ 533 /
|
||||
DATA ( MM( 58, J ), J = 1, 4 ) / 189, 2096, 2159,
|
||||
$ 2841 /
|
||||
DATA ( MM( 59, J ), J = 1, 4 ) / 287, 1761, 2318,
|
||||
$ 4077 /
|
||||
DATA ( MM( 60, J ), J = 1, 4 ) / 2045, 2810, 2091,
|
||||
$ 721 /
|
||||
DATA ( MM( 61, J ), J = 1, 4 ) / 1227, 566, 3443,
|
||||
$ 2821 /
|
||||
DATA ( MM( 62, J ), J = 1, 4 ) / 2838, 442, 1510,
|
||||
$ 2249 /
|
||||
DATA ( MM( 63, J ), J = 1, 4 ) / 209, 41, 449,
|
||||
$ 2397 /
|
||||
DATA ( MM( 64, J ), J = 1, 4 ) / 2770, 1238, 1956,
|
||||
$ 2817 /
|
||||
DATA ( MM( 65, J ), J = 1, 4 ) / 3654, 1086, 2201,
|
||||
$ 245 /
|
||||
DATA ( MM( 66, J ), J = 1, 4 ) / 3993, 603, 3137,
|
||||
$ 1913 /
|
||||
DATA ( MM( 67, J ), J = 1, 4 ) / 192, 840, 3399,
|
||||
$ 1997 /
|
||||
DATA ( MM( 68, J ), J = 1, 4 ) / 2253, 3168, 1321,
|
||||
$ 3121 /
|
||||
DATA ( MM( 69, J ), J = 1, 4 ) / 3491, 1499, 2271,
|
||||
$ 997 /
|
||||
DATA ( MM( 70, J ), J = 1, 4 ) / 2889, 1084, 3667,
|
||||
$ 1833 /
|
||||
DATA ( MM( 71, J ), J = 1, 4 ) / 2857, 3438, 2703,
|
||||
$ 2877 /
|
||||
DATA ( MM( 72, J ), J = 1, 4 ) / 2094, 2408, 629,
|
||||
$ 1633 /
|
||||
DATA ( MM( 73, J ), J = 1, 4 ) / 1818, 1589, 2365,
|
||||
$ 981 /
|
||||
DATA ( MM( 74, J ), J = 1, 4 ) / 688, 2391, 2431,
|
||||
$ 2009 /
|
||||
DATA ( MM( 75, J ), J = 1, 4 ) / 1407, 288, 1113,
|
||||
$ 941 /
|
||||
DATA ( MM( 76, J ), J = 1, 4 ) / 634, 26, 3922,
|
||||
$ 2449 /
|
||||
DATA ( MM( 77, J ), J = 1, 4 ) / 3231, 512, 2554,
|
||||
$ 197 /
|
||||
DATA ( MM( 78, J ), J = 1, 4 ) / 815, 1456, 184,
|
||||
$ 2441 /
|
||||
DATA ( MM( 79, J ), J = 1, 4 ) / 3524, 171, 2099,
|
||||
$ 285 /
|
||||
DATA ( MM( 80, J ), J = 1, 4 ) / 1914, 1677, 3228,
|
||||
$ 1473 /
|
||||
DATA ( MM( 81, J ), J = 1, 4 ) / 516, 2657, 4012,
|
||||
$ 2741 /
|
||||
DATA ( MM( 82, J ), J = 1, 4 ) / 164, 2270, 1921,
|
||||
$ 3129 /
|
||||
DATA ( MM( 83, J ), J = 1, 4 ) / 303, 2587, 3452,
|
||||
$ 909 /
|
||||
DATA ( MM( 84, J ), J = 1, 4 ) / 2144, 2961, 3901,
|
||||
$ 2801 /
|
||||
DATA ( MM( 85, J ), J = 1, 4 ) / 3480, 1970, 572,
|
||||
$ 421 /
|
||||
DATA ( MM( 86, J ), J = 1, 4 ) / 119, 1817, 3309,
|
||||
$ 4073 /
|
||||
DATA ( MM( 87, J ), J = 1, 4 ) / 3357, 676, 3171,
|
||||
$ 2813 /
|
||||
DATA ( MM( 88, J ), J = 1, 4 ) / 837, 1410, 817,
|
||||
$ 2337 /
|
||||
DATA ( MM( 89, J ), J = 1, 4 ) / 2826, 3723, 3039,
|
||||
$ 1429 /
|
||||
DATA ( MM( 90, J ), J = 1, 4 ) / 2332, 2803, 1696,
|
||||
$ 1177 /
|
||||
DATA ( MM( 91, J ), J = 1, 4 ) / 2089, 3185, 1256,
|
||||
$ 1901 /
|
||||
DATA ( MM( 92, J ), J = 1, 4 ) / 3780, 184, 3715,
|
||||
$ 81 /
|
||||
DATA ( MM( 93, J ), J = 1, 4 ) / 1700, 663, 2077,
|
||||
$ 1669 /
|
||||
DATA ( MM( 94, J ), J = 1, 4 ) / 3712, 499, 3019,
|
||||
$ 2633 /
|
||||
DATA ( MM( 95, J ), J = 1, 4 ) / 150, 3784, 1497,
|
||||
$ 2269 /
|
||||
DATA ( MM( 96, J ), J = 1, 4 ) / 2000, 1631, 1101,
|
||||
$ 129 /
|
||||
DATA ( MM( 97, J ), J = 1, 4 ) / 3375, 1925, 717,
|
||||
$ 1141 /
|
||||
DATA ( MM( 98, J ), J = 1, 4 ) / 1621, 3912, 51,
|
||||
$ 249 /
|
||||
DATA ( MM( 99, J ), J = 1, 4 ) / 3090, 1398, 981,
|
||||
$ 3917 /
|
||||
DATA ( MM( 100, J ), J = 1, 4 ) / 3765, 1349, 1978,
|
||||
$ 2481 /
|
||||
DATA ( MM( 101, J ), J = 1, 4 ) / 1149, 1441, 1813,
|
||||
$ 3941 /
|
||||
DATA ( MM( 102, J ), J = 1, 4 ) / 3146, 2224, 3881,
|
||||
$ 2217 /
|
||||
DATA ( MM( 103, J ), J = 1, 4 ) / 33, 2411, 76,
|
||||
$ 2749 /
|
||||
DATA ( MM( 104, J ), J = 1, 4 ) / 3082, 1907, 3846,
|
||||
$ 3041 /
|
||||
DATA ( MM( 105, J ), J = 1, 4 ) / 2741, 3192, 3694,
|
||||
$ 1877 /
|
||||
DATA ( MM( 106, J ), J = 1, 4 ) / 359, 2786, 1682,
|
||||
$ 345 /
|
||||
DATA ( MM( 107, J ), J = 1, 4 ) / 3316, 382, 124,
|
||||
$ 2861 /
|
||||
DATA ( MM( 108, J ), J = 1, 4 ) / 1749, 37, 1660,
|
||||
$ 1809 /
|
||||
DATA ( MM( 109, J ), J = 1, 4 ) / 185, 759, 3997,
|
||||
$ 3141 /
|
||||
DATA ( MM( 110, J ), J = 1, 4 ) / 2784, 2948, 479,
|
||||
$ 2825 /
|
||||
DATA ( MM( 111, J ), J = 1, 4 ) / 2202, 1862, 1141,
|
||||
$ 157 /
|
||||
DATA ( MM( 112, J ), J = 1, 4 ) / 2199, 3802, 886,
|
||||
$ 2881 /
|
||||
DATA ( MM( 113, J ), J = 1, 4 ) / 1364, 2423, 3514,
|
||||
$ 3637 /
|
||||
DATA ( MM( 114, J ), J = 1, 4 ) / 1244, 2051, 1301,
|
||||
$ 1465 /
|
||||
DATA ( MM( 115, J ), J = 1, 4 ) / 2020, 2295, 3604,
|
||||
$ 2829 /
|
||||
DATA ( MM( 116, J ), J = 1, 4 ) / 3160, 1332, 1888,
|
||||
$ 2161 /
|
||||
DATA ( MM( 117, J ), J = 1, 4 ) / 2785, 1832, 1836,
|
||||
$ 3365 /
|
||||
DATA ( MM( 118, J ), J = 1, 4 ) / 2772, 2405, 1990,
|
||||
$ 361 /
|
||||
DATA ( MM( 119, J ), J = 1, 4 ) / 1217, 3638, 2058,
|
||||
$ 2685 /
|
||||
DATA ( MM( 120, J ), J = 1, 4 ) / 1822, 3661, 692,
|
||||
$ 3745 /
|
||||
DATA ( MM( 121, J ), J = 1, 4 ) / 1245, 327, 1194,
|
||||
$ 2325 /
|
||||
DATA ( MM( 122, J ), J = 1, 4 ) / 2252, 3660, 20,
|
||||
$ 3609 /
|
||||
DATA ( MM( 123, J ), J = 1, 4 ) / 3904, 716, 3285,
|
||||
$ 3821 /
|
||||
DATA ( MM( 124, J ), J = 1, 4 ) / 2774, 1842, 2046,
|
||||
$ 3537 /
|
||||
DATA ( MM( 125, J ), J = 1, 4 ) / 997, 3987, 2107,
|
||||
$ 517 /
|
||||
DATA ( MM( 126, J ), J = 1, 4 ) / 2573, 1368, 3508,
|
||||
$ 3017 /
|
||||
DATA ( MM( 127, J ), J = 1, 4 ) / 1148, 1848, 3525,
|
||||
$ 2141 /
|
||||
DATA ( MM( 128, J ), J = 1, 4 ) / 545, 2366, 3801,
|
||||
$ 1537 /
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
I1 = ISEED( 1 )
|
||||
I2 = ISEED( 2 )
|
||||
I3 = ISEED( 3 )
|
||||
I4 = ISEED( 4 )
|
||||
*
|
||||
DO 10 I = 1, MIN( N, LV )
|
||||
*
|
||||
20 CONTINUE
|
||||
*
|
||||
* Multiply the seed by i-th power of the multiplier modulo 2**48
|
||||
*
|
||||
IT4 = I4*MM( I, 4 )
|
||||
IT3 = IT4 / IPW2
|
||||
IT4 = IT4 - IPW2*IT3
|
||||
IT3 = IT3 + I3*MM( I, 4 ) + I4*MM( I, 3 )
|
||||
IT2 = IT3 / IPW2
|
||||
IT3 = IT3 - IPW2*IT2
|
||||
IT2 = IT2 + I2*MM( I, 4 ) + I3*MM( I, 3 ) + I4*MM( I, 2 )
|
||||
IT1 = IT2 / IPW2
|
||||
IT2 = IT2 - IPW2*IT1
|
||||
IT1 = IT1 + I1*MM( I, 4 ) + I2*MM( I, 3 ) + I3*MM( I, 2 ) +
|
||||
$ I4*MM( I, 1 )
|
||||
IT1 = MOD( IT1, IPW2 )
|
||||
*
|
||||
* Convert 48-bit integer to a real number in the interval (0,1)
|
||||
*
|
||||
X( I ) = R*( REAL( IT1 )+R*( REAL( IT2 )+R*( REAL( IT3 )+R*
|
||||
$ REAL( IT4 ) ) ) )
|
||||
*
|
||||
IF (X( I ).EQ.1.0) THEN
|
||||
* If a real number has n bits of precision, and the first
|
||||
* n bits of the 48-bit integer above happen to be all 1 (which
|
||||
* will occur about once every 2**n calls), then X( I ) will
|
||||
* be rounded to exactly 1.0. In IEEE single precision arithmetic,
|
||||
* this will happen relatively often since n = 24.
|
||||
* Since X( I ) is not supposed to return exactly 0.0 or 1.0,
|
||||
* the statistically correct thing to do in this situation is
|
||||
* simply to iterate again.
|
||||
* N.B. the case X( I ) = 0.0 should not be possible.
|
||||
I1 = I1 + 2
|
||||
I2 = I2 + 2
|
||||
I3 = I3 + 2
|
||||
I4 = I4 + 2
|
||||
GOTO 20
|
||||
END IF
|
||||
*
|
||||
10 CONTINUE
|
||||
*
|
||||
* Return final value of seed
|
||||
*
|
||||
ISEED( 1 ) = IT1
|
||||
ISEED( 2 ) = IT2
|
||||
ISEED( 3 ) = IT3
|
||||
ISEED( 4 ) = IT4
|
||||
RETURN
|
||||
*
|
||||
* End of SLARUV
|
||||
*
|
||||
END
|
||||
@ -308,7 +308,7 @@ c ----- print *, itst,iq,t(k),t1a,t1b,cc0,cc1,gvel
|
||||
if(iverb(ifunc).eq.0)then
|
||||
iverb(ifunc) = 1
|
||||
write(LOT,*)'improper initial value in disper - no zero found'
|
||||
write(*,*)'WARNING:improper initial value in disper - no zero found'
|
||||
!write(*,*)'WARNING:improper initial value in disper - no zero found'
|
||||
write(LOT,*)'in fundamental mode '
|
||||
write(LOT,*)'This may be due to low velocity zone '
|
||||
write(LOT,*)'causing reverse phase velocity dispersion, '
|
||||
|
||||
Loading…
Reference in New Issue
Block a user