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

@@ -1,15 +1,15 @@
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c INPUT PARAMETERS
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-surfdataTB.dat			 c: data file 
+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/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)