Update manual again

Fix some bugs in jointly invert both phase and group velocity

bin/DSurfTomo

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

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

@@ -1,21 +1,21 @@
 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 
-4.0  0.0                         c: weight damp
-3                                c: nsublayer (numbers of sublayers for each grid interval:grid --> layer)
-0.5 2.8                          c: minimum velocity, maximum velocity
-10                               c: maxiter (iteration number)
-0.2                              c: sparsity fraction
+20                               c: max(sources, receivers)
+4.0  0.1                         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
+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
 0                                c: kmaxRg
 0                                c: kmaxLc
-0                                c: kmaxLg 
+0                                c: kmaxLg
 0                                c: synthetic flag(0:real data,1:synthetic)
-0.02                             c: noiselevel 
-2.5				 c: threshold
+0.02                             c: noiselevel
+0.5                              c: threshold

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 &

scripts/GenerateIniMOD.py

@@ -6,20 +6,19 @@ import numpy as np
 
 #parameters need to be changed
 #start
-nx=75
-ny=96
-nz=17
+nx=18
+ny=18
+nz=8
 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])
-vel1=np.loadtxt('mod.1d')
+dep1=np.array([0,0.2,0.4,0.6,0.8,1.1,1.4,1.8,2.5])
 #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])

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)
-y=range(1,ny+1)
+nmax = np.max([nx,ny])
+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/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()

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)
@@ -1088,7 +1118,7 @@ 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

src/main.f90

@@ -121,7 +121,7 @@ program SurfTomo
 
   ! OUTPUT PROGRAM INFOMATION            
   write(*,*)
-  write(*,*),'                         S U R F  T O M O'
+  write(*,*),'                             DSurfTomo (v1.3)'
   write(*,*),'PLEASE contact Hongjain Fang &
     (fanghj@mail.ustc.edu.cn) if you find any bug'
   write(*,*)
@@ -181,9 +181,9 @@ program SurfTomo
     if (checkstat > 0) stop 'error allocating RP'
     read(10,*)(tRc(i),i=1,kmaxRc)
     write(*,*)'Rayleigh wave phase velocity used,periods:(s)'
-    write(*,'(50f6.2)')(tRc(i),i=1,kmaxRc)
+    write(*,'(50f7.1)')(tRc(i),i=1,kmaxRc)
     write(66,*)'Rayleigh wave phase velocity used,periods:(s)'
-    write(66,'(50f6.2)')(tRc(i),i=1,kmaxRc)
+    write(66,'(50f7.1)')(tRc(i),i=1,kmaxRc)
   endif
   read(10,*)kmaxRg
   if(kmaxRg.gt.0)then
@@ -191,9 +191,9 @@ program SurfTomo
     if (checkstat > 0) stop 'error allocating RP'
     read(10,*)(tRg(i),i=1,kmaxRg)
     write(*,*)'Rayleigh wave group velocity used,periods:(s)'
-    write(*,'(50f6.2)')(tRg(i),i=1,kmaxRg)
+    write(*,'(50f7.1)')(tRg(i),i=1,kmaxRg)
     write(66,*)'Rayleigh wave group velocity used,periods:(s)'
-    write(66,'(50f6.2)')(tRg(i),i=1,kmaxRg)
+    write(66,'(50f7.1)')(tRg(i),i=1,kmaxRg)
   endif
   read(10,*)kmaxLc
   if(kmaxLc.gt.0)then
@@ -201,9 +201,9 @@ program SurfTomo
     if (checkstat > 0) stop 'error allocating RP'
     read(10,*)(tLc(i),i=1,kmaxLc)
     write(*,*)'Love wave phase velocity used,periods:(s)'
-    write(*,'(50f6.2)')(tLc(i),i=1,kmaxLc)
+    write(*,'(50f7.1)')(tLc(i),i=1,kmaxLc)
     write(66,*)'Love wave phase velocity used,periods:(s)'
-    write(66,'(50f6.2)')(tLc(i),i=1,kmaxLc)
+    write(66,'(50f7.1)')(tLc(i),i=1,kmaxLc)
   endif
   read(10,*)kmaxLg
   if(kmaxLg.gt.0)then
@@ -211,9 +211,9 @@ program SurfTomo
     if (checkstat > 0) stop 'error allocating RP'
     read(10,*)(tLg(i),i=1,kmaxLg)
     write(*,*)'Love wave group velocity used,periods:(s)'
-    write(*,'(50f6.2)')(tLg(i),i=1,kmaxLg)
+    write(*,'(50f7.1)')(tLg(i),i=1,kmaxLg)
     write(66,*)'Love wave group velocity used,periods:(s)'
-    write(66,'(50f6.2)')(tLg(i),i=1,kmaxLg)
+    write(66,'(50f7.1)')(tLg(i),i=1,kmaxLg)
   endif
   read(10,*)ifsyn
   read(10,*)noiselevel
@@ -333,7 +333,7 @@ program SurfTomo
   enddo
   close(10)
   write(*,*) 'grid points in depth direction:(km)'
-  write(*,'(50f6.2)') depz
+  write(*,'(50f7.1)') depz
 
 
 
@@ -382,12 +382,36 @@ program SurfTomo
       cbst(i) = obst(i) - dsyn(i)
     enddo
 
-    threshold = threshold0+(maxiter/2-iter)/3*0.5
+    ! write out rw, iw and cbst for testing
+    !open(44,file='iw.dat')
+    !write(44,*) nar
+    !do i = 2,nar+1
+    !write(44,*) iw(i)
+    !enddo
+    !do i = 1,nar
+    !write(44,*) col(i)
+    !enddo
+    !close(44)
+
+    !open(44,file='rw.dat')
+    !do i = 1,nar
+    !write(44,*) rw(i)
+    !enddo
+    !close(44)
+
+    !open(44,file='residal.dat')
+    !do i = 1,dall
+    !write(44,*) cbst(i)
+    !enddo
+    !close(44)
+
+    !threshold = threshold0+(maxiter/2-iter)/3*0.5
     do i = 1,dall
-      datweight(i) = 1.0
-      if(abs(cbst(i)) > threshold) then
-        datweight(i) = exp(-(abs(cbst(i))-threshold))
-      endif
+!      datweight(i) = 1.0
+!      if(abs(cbst(i)) > threshold) then
+!        datweight(i) = exp(-(abs(cbst(i))-threshold))
+!      endif
+      datweight(i) = 0.01+1.0/(1+0.05*exp(cbst(i)**2*threshold0))
       cbst(i) = cbst(i)*datweight(i)
     enddo
 
@@ -429,7 +453,8 @@ program SurfTomo
 
 
     ! ADDING REGULARIZATION TERM
-    weight=dnrm2(dall,cbst,1)**2/dall*weight0
+    !weight=dnrm2(dall,cbst,1)**2/dall*weight0
+    weight=weight0
     nar_tmp=nar
     nars=0
 

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