diff --git a/example/.gmtdefaults4 b/example/.gmtdefaults4
index dd2074c..3d0aae9 100644
--- a/example/.gmtdefaults4
+++ b/example/.gmtdefaults4
@@ -1,5 +1,5 @@
 #
-#	GMT-SYSTEM 4.5.9 [64-bit] Defaults file
+#	GMT-SYSTEM 4.5.15 [64-bit] Defaults file
 #
 #-------- Plot Media Parameters -------------
 PAGE_COLOR		= white
diff --git a/example/DSurfTomo.in b/example/DSurfTomo.in
index 67f1e4c..7cafb5a 100644
--- a/example/DSurfTomo.in
+++ b/example/DSurfTomo.in
@@ -19,5 +19,3 @@ surfdataTB.dat			 c: data file
 0                                c: synthetic flag(0:real data,1:synthetic)
 0.02                             c: noiselevel 
 2.5				 c: threshold
-1				 c: modest (1: estimate model variation, 0: no estimation)
-30				 c: number of random models
diff --git a/src/CalSurfG.f90 b/src/CalSurfG.f90
index 85b0de6..92e03d9 100644
--- a/src/CalSurfG.f90
+++ b/src/CalSurfG.f90
@@ -1145,7 +1145,8 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
       do ig = 1,igroup
       if (ig ==2 .and. wavetype(srcnum,knumi) == 2) then
         velf(1:nx*ny) = pvRc(1:nx*ny,periods(srcnum,knumi))
-      else
+      endif
+      if (ig ==2 .and. wavetype(srcnum,knumi) == 1) then
         velf(1:nx*ny) = pvLc(1:nx*ny,periods(srcnum,knumi))
       endif
       call gridder(velf)
diff --git a/src/DSurfTomo b/src/DSurfTomo
index 41d8d07..666f79c 100755
Binary files a/src/DSurfTomo and b/src/DSurfTomo differ
diff --git a/src/main.f90 b/src/main.f90
index fc0f8e7..eb65be5 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -218,8 +218,8 @@ program SurfTomo
   read(10,*)ifsyn
   read(10,*)noiselevel
   read(10,*) threshold0
-  read(10,*) modest
-  read(10,*) numrand
+!  read(10,*) modest
+!  read(10,*) numrand
   close(10)
   nrc=nsrc
   kmax=kmaxRc+kmaxRg+kmaxLc+kmaxLg
@@ -607,69 +607,69 @@ program SurfTomo
   close(nout) !close lsmr.txt
   close(66) !close surf_tomo.log
 
-! USE RANDOM MODEL TO OBTAIN THE MODEL VARIATION
-  !modest = 1
-  if (modest ==1) then
-
-    write(*,*) 'model variation estimation begin...'
-    do iter = 1,numrand
-      call init_random_seed()
-      vsftrue=vsf
-      DO K=1,NZ-1
-        DO J=2,NY-1
-          DO I=2,NX-1
-            idx = (k-1)*(ny-2)*(nx-2)+(j-2)*(nx-2)+i-1
-            dv(idx) = 0.1/EXP(2*NORM(idx)/maxnorm)*gaussian()
-            VSFTRUE(I,J,K) = VSF(I,J,K)+dv(idx)
-          ENDDO
-        ENDDO
-      ENDDO
-      write(*,*),'maximum and minimum velocity variation',maxval(dv),minval(dv)
-
-      call synthetic(nx,ny,nz,maxvp,vsftrue,dsyn,&
-        goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
-        tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
-        scxf,sczf,rcxf,rczf,nrc1,nsrc1,kmax,&
-        nsrc,nrc,0.0)
-
-      do i = 1,dall
-        cbst(i) = obst(i) - dsyn(i)
-      enddo
-
-      write(*,*), dnrm2(dall,cbst,1)/sqrt(real(dall)), 1.05*std_devs
-      if (dnrm2(dall,cbst,1)/sqrt(real(dall)) < 1.05*std_devs) then
-        counte = counte + 1
-        modstat(counte,:) = dv
-      endif
-
-    enddo ! iteration for random models
-
-    write(*,*),'number of of models satisfy requirements',counte
-    modsig = 1.0
-    if (counte>0) then
-      do i=1,maxvp
-        !statis
-        !mean = sum(cbst(1:dall))/dall
-        !std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
-        mean = sum(modstat(1:counte,i))/counte
-        stdvs = sqrt(sum(modstat(1:counte,i)**2)/counte-mean**2) 
-        modsig(i) = stdvs
-      enddo
-    endif
-
-    write(*,*),'write model variation to "model_variation.dat"'
-    open (64,file='model_variation.dat')
-    do k=1,nz-1
-      do j=1,ny-2
-        do i=1,nx-2
-          idx = (k-1)*(ny-2)*(nx-2)+(j-1)*(nx-2)+i
-          write(64,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),modsig(idx)
-        enddo
-      enddo
-    enddo
-    close(64)
-    write(*,*) 'finishing model variation estimation'
-  endif
+!! USE RANDOM MODEL TO OBTAIN THE MODEL VARIATION
+!  !modest = 1
+!  if (modest ==1) then
+!
+!    write(*,*) 'model variation estimation begin...'
+!    do iter = 1,numrand
+!      call init_random_seed()
+!      vsftrue=vsf
+!      DO K=1,NZ-1
+!        DO J=2,NY-1
+!          DO I=2,NX-1
+!            idx = (k-1)*(ny-2)*(nx-2)+(j-2)*(nx-2)+i-1
+!            dv(idx) = 0.1/EXP(2*NORM(idx)/maxnorm)*gaussian()
+!            VSFTRUE(I,J,K) = VSF(I,J,K)+dv(idx)
+!          ENDDO
+!        ENDDO
+!      ENDDO
+!      write(*,*),'maximum and minimum velocity variation',maxval(dv),minval(dv)
+!
+!      call synthetic(nx,ny,nz,maxvp,vsftrue,dsyn,&
+!        goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
+!        tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
+!        scxf,sczf,rcxf,rczf,nrc1,nsrc1,kmax,&
+!        nsrc,nrc,0.0)
+!
+!      do i = 1,dall
+!        cbst(i) = obst(i) - dsyn(i)
+!      enddo
+!
+!      write(*,*), dnrm2(dall,cbst,1)/sqrt(real(dall)), 1.05*std_devs
+!      if (dnrm2(dall,cbst,1)/sqrt(real(dall)) < 1.05*std_devs) then
+!        counte = counte + 1
+!        modstat(counte,:) = dv
+!      endif
+!
+!    enddo ! iteration for random models
+!
+!    write(*,*),'number of of models satisfy requirements',counte
+!    modsig = 1.0
+!    if (counte>0) then
+!      do i=1,maxvp
+!        !statis
+!        !mean = sum(cbst(1:dall))/dall
+!        !std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
+!        mean = sum(modstat(1:counte,i))/counte
+!        stdvs = sqrt(sum(modstat(1:counte,i)**2)/counte-mean**2) 
+!        modsig(i) = stdvs
+!      enddo
+!    endif
+!
+!    write(*,*),'write model variation to "model_variation.dat"'
+!    open (64,file='model_variation.dat')
+!    do k=1,nz-1
+!      do j=1,ny-2
+!        do i=1,nx-2
+!          idx = (k-1)*(ny-2)*(nx-2)+(j-1)*(nx-2)+i
+!          write(64,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),modsig(idx)
+!        enddo
+!      enddo
+!    enddo
+!    close(64)
+!    write(*,*) 'finishing model variation estimation'
+!  endif