correcting the readme

before vorotomo

.gitignore

@@ -0,0 +1,7 @@
+*.o
+*.mod
+Applications
+example_*
+.gitignore
+temp
+DSurfTomo

README.md

@@ -1,15 +1,8 @@
-#DSurfTomo
+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.
-DSurfTomo is the surface wave inversion program which can directly invert surface wave 
+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. 
-dispersion data to 3D shear wave speed without the intermediate step of constructing 
+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.
-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 
+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. 
-avoids the assumption of great-circle propagation that is used in most surface wave 
+
-tomographic studies, but which is not appropriate in complex media.
+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 
-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.
 
-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.

bin/DSurfTomo

@@ -0,0 +1 @@
+../src/DSurfTomo

configure

@@ -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 ('--------------------------------------')
-print 'Finishing DSurfTomo compiling'
+print ('DSurfTomo compiling Finished')
-print '--------------------------------------'
+print ('--------------------------------------')

example/.gmtcommands4

@@ -1,10 +0,0 @@
-# GMT common arguments shelf
--Ba0.2f0.1:S velocity (km/s):
--JM2i
--R121.35/121.605/24.975/25.2
--x2.3i
--Y5i
--y-2.7i
--jM2i
-EOF
-

example/.gmtdefaults4

@@ -1,107 +0,0 @@
-#
-#	GMT-SYSTEM 4.5.15 [64-bit] Defaults file
-#
-#-------- Plot Media Parameters -------------
-PAGE_COLOR		= white
-PAGE_ORIENTATION	= landscape
-PAPER_MEDIA		= a4
-#-------- Basemap Annotation Parameters ------
-ANNOT_MIN_ANGLE		= 20
-ANNOT_MIN_SPACING	= 0
-ANNOT_FONT_PRIMARY	= Helvetica
-ANNOT_FONT_SIZE_PRIMARY	= 6p
-ANNOT_OFFSET_PRIMARY	= 0.2c
-ANNOT_FONT_SECONDARY	= Helvetica
-ANNOT_FONT_SIZE_SECONDARY	= 16p
-ANNOT_OFFSET_SECONDARY	= 0.2c
-DEGREE_SYMBOL		= ring
-HEADER_FONT		= Helvetica
-HEADER_FONT_SIZE	= 36p
-HEADER_OFFSET		= 0.5c
-LABEL_FONT		= Helvetica
-LABEL_FONT_SIZE		= 5p
-LABEL_OFFSET		= 0.1c
-OBLIQUE_ANNOTATION	= 1
-PLOT_CLOCK_FORMAT	= hh:mm:ss
-PLOT_DATE_FORMAT	= yyyy-mm-dd
-PLOT_DEGREE_FORMAT	= ddd:mm:ss
-Y_AXIS_TYPE		= hor_text
-#-------- Basemap Layout Parameters ---------
-BASEMAP_AXES		= WESN
-BASEMAP_FRAME_RGB	= black
-BASEMAP_TYPE		= fancy
-FRAME_PEN		= 1.25p
-FRAME_WIDTH		= 0.1c
-GRID_CROSS_SIZE_PRIMARY	= 0c
-GRID_PEN_PRIMARY	= 0.25p
-GRID_CROSS_SIZE_SECONDARY	= 0c
-GRID_PEN_SECONDARY	= 0.5p
-MAP_SCALE_HEIGHT	= 0.2c
-POLAR_CAP		= 85/90
-TICK_LENGTH		= 0.1c
-TICK_PEN		= 0.3p
-X_AXIS_LENGTH		= 25c
-Y_AXIS_LENGTH		= 15c
-X_ORIGIN		= 2.5c
-Y_ORIGIN		= 2.5c
-UNIX_TIME		= FALSE
-UNIX_TIME_POS		= BL/-2c/-2c
-UNIX_TIME_FORMAT	= %Y %b %d %H:%M:%S
-#-------- Color System Parameters -----------
-COLOR_BACKGROUND	= black
-COLOR_FOREGROUND	= white
-COLOR_NAN		= 128
-COLOR_IMAGE		= adobe
-COLOR_MODEL		= rgb
-HSV_MIN_SATURATION	= 1
-HSV_MAX_SATURATION	= 0.1
-HSV_MIN_VALUE		= 0.3
-HSV_MAX_VALUE		= 1
-#-------- PostScript Parameters -------------
-CHAR_ENCODING		= ISOLatin1+
-DOTS_PR_INCH		= 300
-GLOBAL_X_SCALE		= 1
-GLOBAL_Y_SCALE		= 1
-N_COPIES		= 1
-PS_COLOR		= rgb
-PS_IMAGE_COMPRESS	= lzw
-PS_IMAGE_FORMAT		= ascii
-PS_LINE_CAP		= butt
-PS_LINE_JOIN		= miter
-PS_MITER_LIMIT		= 35
-PS_VERBOSE		= FALSE
-TRANSPARENCY		= 0
-#-------- I/O Format Parameters -------------
-D_FORMAT		= %.12lg
-FIELD_DELIMITER		= tab
-GRIDFILE_FORMAT		= nf
-GRIDFILE_SHORTHAND	= FALSE
-INPUT_CLOCK_FORMAT	= hh:mm:ss
-INPUT_DATE_FORMAT	= yyyy-mm-dd
-IO_HEADER		= FALSE
-N_HEADER_RECS		= 1
-NAN_RECORDS		= pass
-OUTPUT_CLOCK_FORMAT	= hh:mm:ss
-OUTPUT_DATE_FORMAT	= yyyy-mm-dd
-OUTPUT_DEGREE_FORMAT	= D
-XY_TOGGLE		= FALSE
-#-------- Projection Parameters -------------
-ELLIPSOID		= WGS-84
-MAP_SCALE_FACTOR	= default
-MEASURE_UNIT		= cm
-#-------- Calendar/Time Parameters ----------
-TIME_FORMAT_PRIMARY	= full
-TIME_FORMAT_SECONDARY	= full
-TIME_EPOCH		= 2000-01-01T12:00:00
-TIME_IS_INTERVAL	= OFF
-TIME_INTERVAL_FRACTION	= 0.5
-TIME_LANGUAGE		= us
-TIME_UNIT		= d
-TIME_WEEK_START		= Sunday
-Y2K_OFFSET_YEAR		= 1950
-#-------- Miscellaneous Parameters ----------
-HISTORY			= TRUE
-INTERPOLANT		= akima
-LINE_STEP		= 0.025c
-VECTOR_SHAPE		= 0
-VERBOSE			= FALSE

example/plotslicenewVs.gmt

@@ -1,49 +0,0 @@
-#!/bin/csh
-
-#-----------------------------------------------------------
-# 2015-06-08 Hongjian Fang
-# example:csh plotslicenewVs.gmt DSurfTomo.inMeasure.dat 0.2 0.4 0.8 1.1
-#-----------------------------------------------------------
-
-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
-set inp3D  = $1
-set ps   = horizVs.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
-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
-
-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
-
-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
-
-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
-
-rm *cpt *grd
-ps2eps -f $ps
-evince *eps &

example_smoothing_clean/DSurfTomo.in

@@ -5,11 +5,11 @@ surfdataTB.dat			 c: data file
 18 18 9                          c: nx ny nz (grid number in lat lon and depth direction)
 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)
-25                              c: nsrc*maxf 
+20                               c: max(sources, receivers)
-4.0  0.0                         c: weight damp
+4.0  1.0                         c: weight damp
-3                                c: nsublayer (numbers of sublayers for each grid interval:grid --> layer)
+3                                c: sablayers (for computing depth kernel, 2~5)
-0.5 2.8                          c: minimum velocity, maximum velocity
+0.5 2.8                          c: minimum velocity, maximum velocity (a priori information)
-10                               c: maxiter (iteration number)
+10                               c: maximum iteration
 0.2                              c: sparsity fraction 
 26                               c: kmaxRc (followed by periods)
 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0
@@ -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
-2.5				 c: threshold
+3.0                              c: threshold

scripts/GenerateDSurfTomoInputFile.py

@@ -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')

scripts/GenerateIniMOD.py

@@ -6,20 +6,20 @@ import numpy as np
 
 #parameters need to be changed
 #start
-nx=75
+nx=18
-ny=96
+ny=18
-nz=17
 minvel=0.8
 velgrad=0.5
-dep1=1.5+np.array([-1.5, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,10.0,11.0,13.0,16.0,20.0,30.0])
+dep1=np.array([0,0.2,0.4,0.6,0.8,1.1,1.4,1.8,2.5])
-vel1=np.loadtxt('mod.1d')
+#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))
 for k in range(nz):
   for j in range(ny):
     for i in range(nx):
-      mod[k*ny+j,i]= vel1[k]/1.75#minvel+dep1[k]*velgrad
+      mod[k*ny+j,i]= minvel+dep1[k]*velgrad
 with open('MOD','w') as fp:
     for i in range(nz):
         fp.write('%5.1f' % dep1[i])
@@ -30,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]),

scripts/GenerateTrueMOD.py

@@ -12,11 +12,12 @@ nz=9
 minvel=0.9
 velgrad=0.6
 dep1=[0,0.2,0.4,0.6,0.8,1.1,1.4,1.8,2.5]
-anosize=1.0
+anosize=0.5
 amplitude=0.4
 #end
-x=range(1,nx+1)
+nmax = np.max([nx,ny])
-y=range(1,ny+1)
+x=range(1,nmax+1)
+y=range(1,nmax+1)
 z=range(1,nz+1)
 z=np.ones(nz)
 bg=np.zeros((nz,nx,ny))
@@ -24,7 +25,7 @@ cross=np.zeros((nz,ny))
 vs1=np.zeros(nz)
 xy=np.kron(np.sin(anosize*np.array(y)),np.sin(anosize*np.array(x)))
 xyz=np.kron(z,xy)
-pxy=xyz.reshape(nz,nx,ny)
+pxy=xyz.reshape(nz,nmax,nmax)
 for k in range(nz):
   for j in range(ny):
     for i in range(nx):
@@ -35,7 +36,7 @@ for k in range(nz):
     for i in range(nx):
 	mod[(k)*ny+j,i]=bg[k,i,j]+pxy[k,i,j]*amplitude
 k=5
-plt.imshow(mod[k*ny:(k+1)*ny,:],cmap='jet_r',interpolation='bicubic')
+plt.imshow(mod[k*ny:(k+1)*ny,:],cmap='jet_r')
 np.savetxt('MOD.true',mod,fmt='%4.4f')
 plt.colorbar()
 plt.show()

scripts/plotcross.gmt

@@ -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
+

scripts/plotpath.py

@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+
+import numpy as np
+from matplotlib import pylab as plt
+
+#parameters you need to change
+#begin
+numray = 3320 # how many rays do you have (all periods)
+rb = 1 # ray index for some period (first one)
+re = 50 # ray index for some period (last one)
+maxseg = 500 # 
+#end
+
+n = 0
+numseg = np.zeros(numray,)
+raylat = np.zeros((numray,maxseg))
+raylon = np.zeros((numray,maxseg))
+with open('raypath.out') as ray:
+    for line in ray:
+        linseg = line.split()
+        if linseg[0]=='#':
+            n = n+1
+            seg = 0
+            numseg[n-1] = float(linseg[1])
+        else:
+            seg = seg+1
+            raylat[n-1,seg-1] = float(linseg[0])
+            raylon[n-1,seg-1] = float(linseg[1])
+
+for ii in range(rb,re):
+    plt.plot(raylon[ii,0:numseg[ii]],raylat[ii,0:numseg[ii]],'k-')
+
+plt.show()

scripts/plotslice.gmt

@@ -7,42 +7,42 @@
 # csh plotslice.gmt SurfTomo.in.tvMeasure.dat 0.2 0.4 0.8 1.4
 #-----------------------------------------------------------
 
-gmtset ANOT_FONT_SIZE 6
+#gmt set ANOT_FONT_SIZE 6
-gmtset LABEL_FONT_SIZE 6
+#gmt set LABEL_FONT_SIZE 6
-gmtset FRAME_WIDTH=0.1c
+#gmt set FRAME_WIDTH=0.1c
-gmtset LABEL_OFFSET=0.1c
+#gmt set LABEL_OFFSET=0.1c
-gmtset LABEL_FONT_SIZE=5p
+#gmt set LABEL_FONT_SIZE=5p
-gmtset TICK_LENGTH=0.1c
+#gmt set TICK_LENGTH=0.1c
-gmtset TICK_PEN=0.3p
+#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
+gmt makecpt -Cseis -D -T0.5/2.5/0.1   > $cpt #velocity boundary
-set RMAP = `minmax -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
+set RMAP = `gmt gmtinfo -C $inp3D | awk '{print "-R"$1"/"$2"/"$3"/"$4}'`
-psbasemap $RMAP $J -Ba0.1f0.05WseN -P -K -Y5i >$ps
+gmt 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 
+awk '{if($3==depth1) print $1,$2,$4}' depth1=$2 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd 
-grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
+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
+#gmt makecpt -Cseis -D -T0.8/1.8/0.1   > $cpt #velocity boundary
-psbasemap $RMAP $J -Ba0.1f0.05NwsE -P -K -O -X2.3i >>$ps
+gmt 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 
+awk '{if($3==depth2) print $1,$2,$4}' depth2=$3 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd 
-grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
+gmt 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 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
+#gmt makecpt -Cseis -D -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
+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|xyz2grd -R -I0.017/0.015 -Gtmp.grd 
+awk '{if($3==depth2) print $1,$2,$4}' depth2=$4 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd 
-grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
+gmt 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 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
+#gmt makecpt -Cseis -D -T1.3/2.4/0.1   > $cpt #velocity boundary
-psbasemap $RMAP $J -Ba0.1f0.05SwnE -P -K -O -X2.3i >>$ps
+gmt 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 
+awk '{if($3==depth2) print $1,$2,$4}' depth2=$5 $inp3D|gmt xyz2grd -R -I0.017/0.015 -Gtmp.grd 
-grdimage tmp.grd $J $RMAP -BSenW -E100 -C$cpt -O -K >> $ps
+gmt 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 psscale -Cslice.cpt -Ba0.2f0.1:'S velocity (km/s)': -D1.0i/-0.25i/3.00/0.2h  -O  -P >> $ps
 

src/CalSurfG.f90

@@ -92,6 +92,17 @@ subroutine depthkernel(nx,ny,nz,vel,pvRc,sen_vsRc,sen_vpRc,sen_rhoRc, &
           dlncg_dlnvs(nn,i) = (cg2(nn)-cg1(nn))/(dlnVs*vsz(i))
         enddo
 
+        !note here, no integral from z1 to z2 is needed. The grid based nodes 
+        ! have already taken into consideration of it. Layer based need integration.
+        ! Test passed
+        !if (i == 1) then
+        !  dlncg_dlnvs(1:kmaxRc,i) = dlncg_dlnvs(1:kmaxRc,i)*thkm(1)/2.0 
+        !else if (i == mmax) then
+        !  dlncg_dlnvs(1:kmaxRc,i) = dlncg_dlnvs(1:kmaxRc,i)*thkm(mmax-1)/2.0 
+        !else
+        !  dlncg_dlnvs(1:kmaxRc,i) = dlncg_dlnvs(1:kmaxRc,i)*(thkm(i-1)+thkm(i))/2.0 
+        !endif
+
 
         vpm(i) = vpz(i) - 0.5*dlnVp*vpz(i)
         call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,&
@@ -110,6 +121,16 @@ subroutine depthkernel(nx,ny,nz,vel,pvRc,sen_vsRc,sen_vpRc,sen_rhoRc, &
           !                cga = 0.5*(cg1(nn)+cg2(nn))
           dlncg_dlnvp(nn,i) = (cg2(nn)-cg1(nn))/(dlnVp*vpz(i))
         enddo
+
+
+        !if (i == 1) then
+        !  dlncg_dlnvp(1:kmaxRc,i) = dlncg_dlnvp(1:kmaxRc,i)*thkm(1)/2.0 
+        !else if (i == mmax) then
+        !  dlncg_dlnvp(1:kmaxRc,i) = dlncg_dlnvp(1:kmaxRc,i)*thkm(mmax-1)/2.0 
+        !else
+        !  dlncg_dlnvp(1:kmaxRc,i) = dlncg_dlnvp(1:kmaxRc,i)*(thkm(i-1)+thkm(i))/2.0 
+        !endif
+
         rhom(i) = rhoz(i) - 0.5*dlnrho*rhoz(i)
         call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,&
           rmax,rdep,rvp,rvs,rrho,rthk)
@@ -127,6 +148,15 @@ subroutine depthkernel(nx,ny,nz,vel,pvRc,sen_vsRc,sen_vpRc,sen_rhoRc, &
           !                cga = 0.5*(cg1(nn)+cg2(nn))
           dlncg_dlnrho(nn,i) = (cg2(nn)-cg1(nn))/(dlnrho*rhoz(i))
         enddo
+
+       ! if (i == 1) then
+       !   dlncg_dlnrho(1:kmaxRc,i) = dlncg_dlnrho(1:kmaxRc,i)*thkm(1)/2.0 
+       ! else if (i == mmax) then
+       !   dlncg_dlnrho(1:kmaxRc,i) = dlncg_dlnrho(1:kmaxRc,i)*thkm(mmax-1)/2.0 
+       ! else
+       !   dlncg_dlnrho(1:kmaxRc,i) = dlncg_dlnrho(1:kmaxRc,i)*(thkm(i-1)+thkm(i))/2.0 
+       ! endif
+
       enddo
       sen_vsRc((jj-1)*nx+ii,1:kmaxRc,1:mmax)=dlncg_dlnvs(1:kmaxRc,1:mmax)
       sen_vpRc((jj-1)*nx+ii,1:kmaxRc,1:mmax)=dlncg_dlnvp(1:kmaxRc,1:mmax)
@@ -910,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
@@ -997,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
+  gdx=8
-  gdz=5
+  gdz=8
   asgr=1
   sgdl=8
   sgs=8
@@ -1069,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
@@ -1078,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
@@ -1088,14 +1119,14 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
     iwave=1
     igr=0
     call depthkernel(nx,ny,nz,vels,pvLc,sen_vsLc,sen_vpLc, &
-      sen_rhoLc,iwave,igr,kmax,tLc,depz,minthk)
+      sen_rhoLc,iwave,igr,kmaxLc,tLc,depz,minthk)
   endif
 
   if(kmaxLg.gt.0) then
     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)
@@ -1252,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
@@ -1346,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
@@ -1370,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
@@ -1716,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)
@@ -1737,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
@@ -2222,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
+  !if(writepath == 1) then
-    WRITE(40,*)'#',nrp
+  !  WRITE(40,*)'#',nrp
-    DO j=1,nrp
+  !  DO j=1,nrp
-      rayx=(pi/2-rgx(j))*180.0/pi
+  !    rayx=(pi/2-rgx(j))*180.0/pi
-      rayz=rgz(j)*180.0/pi
+  !    rayz=rgz(j)*180.0/pi
-      WRITE(40,*)rayx,rayz
+  !    WRITE(40,*)rayx,rayz
-    ENDDO
+  !  ENDDO
-  endif
+  !endif
   !fang!   ENDIF
   !
   !  Write partial derivatives to output file

src/Makefile

@@ -1,17 +1,19 @@
 CMD = DSurfTomo
 FC = gfortran
-FFLAGS  = -O3 -ffixed-line-length-none -ffloat-store\
+FFLAGS  = -O -ffixed-line-length-none -ffloat-store\
-           -W  -fbounds-check -m64 -mcmodel=medium
+           -fbounds-check -m64 -mcmodel=medium
 F90SRCS = lsmrDataModule.f90 lsmrblasInterface.f90\
           lsmrblas.f90 lsmrModule.f90 delsph.f90\
-	  aprod.f90 gaussian.f90 main.f90
+	      aprod.f90 gaussian.f90 getpercentile.f90
-FSRCS =  surfdisp96.f  
+FSRCS =  surfdisp96.f slarnv.f slaruv.f
-OBJS = $(F90SRCS:%.f90=%.o) $(FSRCS:%.f=%.o) CalSurfG.o 
+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

src/getpercentile.f90

@@ -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
+

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
+        integer,dimension(:),allocatable::iw,col,nrow
-  real,dimension(:),allocatable::dv,norm
+        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
-  nout=36
+        !open(nout,file='lsmr.txt')
-  open(nout,file='lsmr.txt')
 
         ! OUTPUT PROGRAM INFOMATION            
         write(*,*)
-  write(*,*),'                         S U R F  T O M O'
+        write(*,*) '                             DSurfTomo (v1.4)'
-  write(*,*),'PLEASE contact Hongjain Fang &
+        !write(*,*) 'PLEASE contact Hongjain Fang &
-    (fanghj@mail.ustc.edu.cn) if you find any bug'
+        !        (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
@@ -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,*)'                         S U R F  T O M O'
-  write(66,*),'PLEASE contact Hongjain Fang &
+        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(*,'(50f6.2)')(tRc(i),i=1,kmaxRc)
+                write(*,'(50f7.2)')(tRc(i),i=1,kmaxRc)
                 write(66,*)'Rayleigh wave phase velocity used,periods:(s)'
-    write(66,'(50f6.2)')(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(*,'(50f6.2)')(tRg(i),i=1,kmaxRg)
+                write(*,'(50f7.2)')(tRg(i),i=1,kmaxRg)
                 write(66,*)'Rayleigh wave group velocity used,periods:(s)'
-    write(66,'(50f6.2)')(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(*,'(50f6.2)')(tLc(i),i=1,kmaxLc)
+                write(*,'(50f7.2)')(tLc(i),i=1,kmaxLc)
                 write(66,*)'Love wave phase velocity used,periods:(s)'
-    write(66,'(50f6.2)')(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(*,'(50f6.2)')(tLg(i),i=1,kmaxLg)
+                write(*,'(50f7.2)')(tLg(i),i=1,kmaxLg)
                 write(66,*)'Love wave group velocity used,periods:(s)'
-    write(66,'(50f6.2)')(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
@@ -333,13 +320,11 @@ program SurfTomo
         enddo
         close(10)
         write(*,*) 'grid points in depth direction:(km)'
-  write(*,'(50f6.2)') 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,36 +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
 
-    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 (cbst(i)<q25*threshold0 .or. cbst(i)>q75*threshold0) then
-      if(abs(cbst(i)) > threshold) then
+            datweight(i) = 0.0
-        datweight(i) = exp(-(abs(cbst(i))-threshold))
+            cbst(i) = 0
         endif
-      cbst(i) = cbst(i)*datweight(i)
         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))
@@ -407,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
@@ -428,8 +411,7 @@ program SurfTomo
         endif
 
 
-    ! ADDING REGULARIZATION TERM
+        weight=weight0
-    weight=dnrm2(dall,cbst,1)**2/dall*weight0
         nar_tmp=nar
         nars=0
 
@@ -489,10 +471,12 @@ program SurfTomo
         leniw = 2*nar+1
         lenrw = nar
         dv = 0
-    atol = 1e-3
+        atol = 1e-6
-    btol = 1e-3
+        btol = 1e-6
-    conlim = 1200
+        !atol = 1e-3/((dvxd+dvzd)*111.19/2.0*0.1) !1e-2
-    itnlim = 1000
+        !btol = 1e-3/(dvxd*nx*111.19/3.0)!1e-3
+        conlim = 100
+        itnlim = 400
         istop = 0
         anorm = 0.0
         acond = 0.0
@@ -503,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(*,'(i2,a)')iter,'th iteration...'
-    write(*,'(a,f7.3)'),'weight is:',weight
+!        write(*,'(a,f7.3)')'weight is:',weight
-    write(*,'(a,f8.1,a,f8.2,a,f8.3)'),'mean,std_devs and rms of &
+        write(*,'(a,f8.1,a,f8.2,a,f8.3)')' mean,std_devs and rms of &
-      residual: ',mean*1000,'ms ',1000*std_devs,'ms ',&
+                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
+        !do i =1,dall
-    write(66,'(a,f8.1,a,f8.2,a,f8.3)'),'mean,std_devs and rms of &
+        !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
@@ -543,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
@@ -572,8 +547,8 @@ program SurfTomo
 
         ! OUTPUT THE VELOCITY MODEL
 
-  write(*,*),'Program finishes successfully'
+        write(*,*)'Program finishes successfully'
-  write(66,*),'Program finishes successfully'
+        write(66,*)'Program finishes successfully'
 
         if(ifsyn == 1) then
                 open(65,file='Vs_model.real')
@@ -582,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(65,'(5f10.5)') 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(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'
@@ -603,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(40)
-  close(34)
+        !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)
@@ -694,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)
@@ -715,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

src/slarnv.f

@@ -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

src/slaruv.f

@@ -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

src/surfdisp96.f

@@ -56,7 +56,7 @@ c-----
         integer NL, NL2, NLAY
         parameter(NL=200,NLAY=200,NL2=NL+NL)
         integer NP
-        parameter (NP=60)
+        parameter (NP=80)
 
 c-----
 c     LIN - unit for FORTRAN read from terminal
@@ -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, '
@@ -500,7 +500,7 @@ c                 2 - Rayleigh Wave
 c     iflag   I*4 0 - Initialize
 c                 1 - Make model  for Love or Rayleigh Wave
 c-----
-        parameter(NL=200,NP=60)
+        parameter(NL=200,NP=80)
         real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
         integer mmax,llw
 c        common/modl/ d,a,b,rho,rtp,dtp,btp
@@ -563,7 +563,7 @@ c     nev = 0 force interval halving
 c     nev = 1 permit neville iteration if conditions are proper
 c     nev = 2 neville iteration is being used
 c-----
-        parameter (NL=200,NP=60)
+        parameter (NL=200,NP=80)
         implicit double precision (a-h,o-z)
         real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
         dimension x(20),y(20)
@@ -704,7 +704,7 @@ c
         function dltar1(wvno,omega,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi)
 c   find SH dispersion values.
 c
-        parameter (NL=200,NP=60)
+        parameter (NL=200,NP=80)
         implicit double precision (a-h,o-z)
         real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)
         integer llw,mmax
@@ -768,7 +768,7 @@ c
 c     &  a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz)
 c   find P-SV dispersion values.
 c
-        parameter (NL=200,NP=60)
+        parameter (NL=200,NP=80)
         implicit double precision (a-h,o-z)
         dimension e(5),ee(5),ca(5,5)
         real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)