From b657c68d2877fa412db44193cf24ff06fb5e4770 Mon Sep 17 00:00:00 2001
From: Hongjian <fanghj1990@gmail.com>
Date: Sun, 19 Apr 2020 07:52:07 -0400
Subject: [PATCH] -

---
 src/CalSurfG.f90    |  8 +++++---
 src/Makefile        |  2 +-
 src/main.f90        | 10 +++++++---
 src/voronoiproj.f90 | 12 ++++++------
 4 files changed, 19 insertions(+), 13 deletions(-)

diff --git a/src/CalSurfG.f90 b/src/CalSurfG.f90
index ce33aa4..788a0a2 100644
--- a/src/CalSurfG.f90
+++ b/src/CalSurfG.f90
@@ -1098,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
@@ -1109,6 +1110,7 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
     call caldespersion(nx,ny,nz,vels,pvRc, &
         iwave,igr,kmax,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
@@ -1377,9 +1379,9 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, &
         !
         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
+        !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))
         row(1:nparpi)=0.0
diff --git a/src/Makefile b/src/Makefile
index 077ff65..f295801 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,6 +1,6 @@
 CMD = DSurfTomo
 FC = gfortran
-FFLAGS  = -O3 -ffixed-line-length-none -ffloat-store\
+FFLAGS  = -O -ffixed-line-length-none -ffloat-store\
            -fbounds-check -m64 -mcmodel=medium
 F90SRCS = lsmrDataModule.f90 lsmrblasInterface.f90\
           lsmrblas.f90 lsmrModule.f90 delsph.f90\
diff --git a/src/main.f90 b/src/main.f90
index d716830..21be9c4 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -138,7 +138,7 @@ program SurfTomo
         read(10,'(a30)')dummy
         read(10,'(a30)')dummy
         read(10,'(a30)')dummy
-        read(10,*)datafile
+        read(10,*) datafile
         read(10,*) nx,ny,nz
         read(10,*) goxd,gozd
         read(10,*) dvxd,dvzd
@@ -289,6 +289,7 @@ program SurfTomo
         close(87)
         allocate(depz(nz), stat=checkstat)
         maxnar = spfra*dall*nx*ny*nz!sparsity fraction
+        if (maxnar<0) print*, 'number overflow, decrease your sparsefrac'
         maxvp = (nx-2)*(ny-2)*(nz-1)
         allocate(dv(maxvp),dvsub(maxvp),dvstd(maxvp),dvall(maxvp*nrealizations), stat=checkstat)
 !        allocate(dvall(maxvp*nrealizations),stats=checkstat)
@@ -416,6 +417,7 @@ program SurfTomo
             lenrw = nar
             iw(1)=nar
             iw(nar+2:2*nar+1) = col(1:nar)
+            print*,'no. of nonzero:',nar,minval(cbst),maxval(cbst)
             !$omp parallel &
             !$omp default(private) &
             !$omp shared(leniw,lenrw,iw,rw,cbst,goxd,gozd,dvxd,dvzd,depz,maxvp) &
@@ -497,8 +499,10 @@ program SurfTomo
         leniw = 2*nar+1
         lenrw = nar
         dv = 0
-        atol = 1e-4
-        btol = 1e-4
+        !atol = 1e-4
+        !btol = 1e-4
+        atol = 1e-3/((dvxd+dvzd)*111.19/2.0*0.1) !1e-2
+        btol = 1e-3/(dvxd*nx*111.19/3.0)!1e-3
         conlim = 100
         itnlim = 400
         istop = 0
diff --git a/src/voronoiproj.f90 b/src/voronoiproj.f90
index 6ef4a65..0f07e5e 100644
--- a/src/voronoiproj.f90
+++ b/src/voronoiproj.f90
@@ -140,8 +140,8 @@ subroutine voronoiproj(leniw,lenrw,iw,rw,dres,goxd,dvxd,gozd,dvzd,depz,&
 
       conlim = 50
       itnlim = 100
-      atol = 1e-2
-      btol = 1e-3
+      atol = 1e-3/((dvxd+dvzd)*111.19/2.0*0.1) !1e-2
+      btol = 1e-3/(dvxd*nx*111.19/3.0)!1e-3
       istop = 0
       anorm = 0.0
       acond = 0.0
@@ -152,14 +152,14 @@ subroutine voronoiproj(leniw,lenrw,iw,rw,dres,goxd,dvxd,gozd,dvzd,depz,&
       ! using lsmr to solve for the projection coefficients
       !print*, 'LSMR beginning ...'
 
-      nout = -1
-      !nout = 36
-      !open(nout,file='lsmrout_sub.txt')
+      !nout = -1
+      nout = 36
+      open(nout,file='lsmrout_sub.txt')
 
       call LSMR(nd, ncells, leniwgp, lenrwgp,iwgp,rwgp,dres,damp,&
       atol, btol, conlim, itnlim, localSize,nout,&
       xunknown, istop, itn, anorm, acond,rnorm, arnorm, xnorm)
-      !close(nout)
+      close(nout)
       do ii = 1,ncells
         xunknown(ii) = xunknown(ii)/norm(ii)
       enddo