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.

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.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_smoothing_clean/DSurfTomo.in

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

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

@@ -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]),

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

@@ -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
+  gdx=8
-  gdz=5
+  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
+  !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/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

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