change configure

README.md

@@ -1 +0,0 @@
-# SurfTomo15

configure

@@ -0,0 +1,16 @@
+#!/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'
+else:
+ print 'installation beginning'
+ os.mkdir('bin')
+os.chdir('src')
+os.system('make clean')
+os.system('make')
+os.system('cp SurfTomo ../bin')
+print '--------------------------------------'
+print 'surf_tomo install over'
+print '--------------------------------------'

src/Makefile

@@ -0,0 +1,25 @@
+CMD = SurfTomo
+FC = gfortran
+FFLAGS  = -O3 -ffixed-line-length-none -ffloat-store\
+           -W  -fbounds-check -m64 -mcmodel=medium
+F90SRCS = lsmrDataModule.f90 lsmrblasInterface.f90  \
+          lsmrblas.f90 lsmrModule.f90 delsph.f90\
+	  forwardstep.f90 forwardtrans.f90 split.f90 merge1.f90\
+	  invtrans3d.f90 inversetrans.f90 aprod.f90\
+	  haar.f90 waveletD8.f90 inversestep.f90\
+	  gaussian.f90 main.f90
+FSRCS =  surfdisp96.f  
+OBJS = $(F90SRCS:%.f90=%.o) $(FSRCS:%.f=%.o) CalSurfG.o wavelettrans3domp.o
+all:$(CMD)
+$(CMD):$(OBJS)
+	$(FC) -fopenmp $^ -o $@
+CalSurfG.o:CalSurfG.f90
+	$(FC) -fopenmp $(FFLAGS) -c $<  -o $@
+wavelettrans3domp.o: wavelettrans3domp.f90
+	$(FC) -fopenmp $(FFLAGS) -c $<  -o $@
+%.o: %.f90
+	$(FC) $(FFLAGS) -c $(@F:.o=.f90) -o $@
+%.o: %.f
+	$(FC) $(FFLAGS) -c $(@F:.o=.f) -o $@
+clean:
+	rm *.o *.mod $(CMD)

src/aprod.f90

@@ -0,0 +1,60 @@
+!c--- This file is from hypoDD by Felix Waldhauser ---------
+!c-------------------------Modified by Haijiang Zhang-------
+!c Multiply a matrix by a vector
+!c  Version for use with sparse matrix specified by
+!c  output of subroutine sparse for use with LSQR
+
+	subroutine aprod(mode, m, n, x, y, leniw, lenrw, iw, rw)
+
+	implicit none
+
+!c	Parameters:
+	integer	mode		! ==1: Compute  y = y + a*x
+				! 	y is altered without changing x
+				! ==2: Compute  x = x + a(transpose)*y
+				!	x is altered without changing y
+	integer	m, n		! Row and column dimensions of a
+	real	x(n), y(m)	! Input vectors
+	integer :: leniw
+	integer lenrw
+	integer	iw(leniw)	! Integer work vector containing:
+				! iw[1]  Number of non-zero elements in a
+				! iw[2:iw[1]+1]  Row indices of non-zero elements
+				! iw[iw[1]+2:2*iw[1]+1]  Column indices
+	real	rw(lenrw)	! [1..iw[1]] Non-zero elements of a
+
+!c	Local variables:
+	integer i1
+	integer j1
+	integer k
+	integer kk
+
+!c	set the ranges the indices in vector iw
+
+	kk=iw(1)
+	i1=1
+	j1=kk+1
+
+!c	main iteration loop
+
+	do k = 1,kk
+
+	if (mode.eq.1) then
+
+!c	compute  y = y + a*x
+
+	y(iw(i1+k)) = y(iw(i1+k)) + rw(k)*x(iw(j1+k))
+
+	else
+
+!c	compute  x = x + a(transpose)*y
+
+	x(iw(j1+k)) = x(iw(j1+k)) + rw(k)*y(iw(i1+k))
+
+	endif
+        enddo
+
+!  100	continue
+
+	return
+	end

src/delsph.f90

@@ -0,0 +1,28 @@
+subroutine delsph(flat1,flon1,flat2,flon2,del)
+implicit none
+real,parameter:: R=6371.0
+REAL,parameter:: pi=3.1415926535898
+real flat1,flat2
+real flon1,flon2
+real del
+
+real dlat
+real dlon
+real lat1
+real lat2
+real a
+real c
+
+
+!dlat=(flat2-flat1)*pi/180
+!dlon=(flon2-flon1)*pi/180
+!lat1=flat1*pi/180
+!lat2=flat2*pi/180
+dlat=flat2-flat1
+dlon=flon2-flon1
+lat1=pi/2-flat1
+lat2=pi/2-flat2
+a=sin(dlat/2)*sin(dlat/2)+sin(dlon/2)*sin(dlon/2)*cos(lat1)*cos(lat2)
+c=2*atan2(sqrt(a),sqrt(1-a))
+del=R*c
+end subroutine

src/forwardstep.f90

@@ -0,0 +1,26 @@
+	subroutine forwardstep(vec,n)
+	implicit none
+	integer n
+	real vec(n)
+	
+	integer half, k, j
+	half=n/2!changed to more than one
+	do k=1,half
+	vec(k)=vec(k)+sqrt(3.0)*vec(half+k)
+	end do
+	vec(half+1)=vec(half+1)-sqrt(3.0)/4.0*vec(1)- &
+        (sqrt(3.0)-2)/4.0*vec(half)
+	do k=1,half-1
+	vec(half+k+1)=vec(half+k+1)-sqrt(3.0)/4.0*vec(k+1)- &
+        (sqrt(3.0)-2)/4.0*vec(k)
+	end do
+	do k=1,half-1
+	vec(k)=vec(k)-vec(half+k+1)
+	end do
+	vec(half)=vec(half)-vec(half+1) 
+	
+	do k=1,half
+	vec(k)=(sqrt(3.0)-1)/sqrt(2.0)*vec(k)
+	vec(half+k)=(sqrt(3.0)+1)/sqrt(2.0)*vec(half+k)
+	end do
+	end subroutine

src/forwardtrans.f90

@@ -0,0 +1,47 @@
+	subroutine forwardtrans(vec,n,mxl,tp)
+	implicit none
+	integer n,mxl,tp
+	real vec(n)
+	integer i,j
+	integer forward
+	i=n
+	if (tp == 1) then
+	forward=0
+	do while(i.ge.n/(2**mxl))
+	call split(vec,i)
+	call predict(vec,i,forward)
+	call update(vec,i,forward)
+        call normalizationf(vec,i)
+	i=i/2
+	enddo
+	endif
+
+	if (tp == 2) then
+        do while(i.ge.n/(2**mxl))
+        call split(vec,i)
+        call forwardstep(vec,i)
+        i=i/2
+        enddo
+	endif
+
+	if (tp == 3) then
+	do while(i.ge.n/(2**mxl))
+	call transformD8(vec,i)
+	i=i/2
+	enddo
+	endif 
+	end subroutine
+
+	subroutine normalizationf(x,n)
+        implicit none
+        real x(*)
+        integer n
+        
+        integer k
+        do k=1,n/2
+        x(k)=x(k)*sqrt(2.0)
+        enddo
+        do k=n/2+1,n
+        x(k)=x(k)*sqrt(2.0)/2
+        enddo
+        end

src/gaussian.f90

@@ -0,0 +1,31 @@
+	real function gaussian()
+	implicit none
+!	real rd
+	
+	real x1,x2,w,y1
+	real y2
+	real n1,n2
+	integer use_last
+	integer ii,jj
+
+	use_last=0
+	y2=0
+	w=2.0
+	if(use_last.ne.0) then
+	  y1=y2
+	  use_last=0
+	else
+	  do while (w.ge.1.0)
+	   call random_number(n1)
+	   call random_number(n2)
+	    x1=2.0*n1-1.0
+	    x2=2.0*n2-1.0
+	    w = x1 * x1 + x2 * x2
+	  enddo
+	w=((-2.0*log(w))/w)**0.5
+	y1=x1*w
+	y2=x2*w
+	use_last=1
+	endif
+	gaussian=y1
+	end function

src/haar.f90

@@ -0,0 +1,49 @@
+  	subroutine predict(vec, N, direction )
+    	implicit none
+	real vec(*)
+	integer N           !must be power of 2
+	integer direction   !0:forward 1:inverse
+	!local variables
+	integer half
+	integer cnt,i,j
+	real predictVal
+    	half = N/2
+    	cnt = 0
+
+	do i=1,half
+        predictVal = vec(i)
+        j = i + half
+        if(direction == 0) then
+	vec(j) = vec(j) - predictVal
+        else if (direction == 1) then
+	vec(j) = vec(j) + predictVal
+        else 
+	print*,"haar::predict: bad direction value"
+        stop
+        endif
+	enddo
+	end subroutine
+	
+    	subroutine update( vec, N, direction )
+    	implicit none
+	real vec(*)
+	integer N           !must be power of 2
+	integer direction   !0:forward 1:inverse
+	!local variables
+	integer half
+	integer cnt,i,j
+	real updateVal
+    	half = N/2
+    	do i=1,half
+        j = i + half
+        updateVal = vec(j) / 2.0
+        if (direction == 0) then
+	vec(i) = vec(i) + updateVal;
+        else if (direction ==1) then
+	vec(i) = vec(i) - updateVal
+        else 
+	print*,"update: bad direction value"
+	stop
+	endif
+	enddo
+	end subroutine

src/inversestep.f90

@@ -0,0 +1,25 @@
+	subroutine inversestep(vec,n)
+	implicit none
+	integer n
+	real vec(n)
+	
+	integer half, k
+	half=int(n/2.0)
+	do k=1,half,1
+	vec(k)=(sqrt(3.0)+1.0)/sqrt(2.0) * vec(k)
+	vec(k+half)=(sqrt(3.0)-1.0)/sqrt(2.0) * vec(k+half)
+	enddo
+	do k=1,half-1,1
+	vec(k)=vec(k)+vec(half+k+1)
+	enddo
+	vec(half)=vec(half)+vec(half+1)
+	vec(half+1)=vec(half+1)+sqrt(3.0)/4.0*vec(1)+ &
+       (sqrt(3.0)-2)/4.0*vec(half)
+	do k=2,half,1
+	vec(half+k)=vec(half+k)+sqrt(3.0)/4.0*vec(k)+ &
+        (sqrt(3.0)-2)/4.0*vec(k-1)
+	enddo
+	do k=1,half,1
+	vec(k)=vec(k)-sqrt(3.0)*vec(half+k)
+	enddo
+	end subroutine

src/inversetrans.f90

@@ -0,0 +1,47 @@
+	subroutine inversetrans(vec,n,mxl,tp)
+	implicit none
+	integer n
+	real vec(n)
+	integer i
+	integer mxl,tp
+	integer inverse
+	i=n/(2**mxl) !n-->n/2 
+	if (tp == 1) then
+	inverse=1
+	do while (i.le.n)
+        call normalizationi(vec,i)
+	call update(vec,i,inverse)
+	call predict(vec,i,inverse)
+	call merge1(vec,i)
+	i=i*2
+	enddo
+	endif
+	if (tp == 2) then
+        do while (i.le.n)
+        call inversestep(vec,i)
+        call merge1(vec,i)
+        i=i*2
+        enddo
+	endif
+	if (tp == 3) then
+	do while (i.le.n)
+	call invTransformD8(vec,i)
+	i=i*2
+	enddo
+	endif
+	end subroutine
+
+
+        subroutine normalizationi(x,i)
+        implicit none
+        real x(*)
+        integer i
+        
+        integer k
+        do k=1,i/2
+        x(k)=x(k)/sqrt(2.0)
+        enddo
+        do k=i/2+1,i
+        x(k)=x(k)*sqrt(2.0)
+        enddo
+        end

src/invtrans3d.f90

@@ -0,0 +1,32 @@
+	subroutine invwavetrans(nx,ny,nz,x,mxl,mxld,HorizonType,VerticalType)
+	implicit none
+	integer mxl,mxld,HorizonType,VerticalType
+	integer nx,ny,nz
+	real x(*)
+	real fxs(nx),fys(ny),fzs(nz)
+	integer k,j,i
+!c 	local variables
+	do k=1,ny
+	do j=1,nx
+	fzs=x(j+(k-1)*nx:j+(k-1)*nx+(nz-1)*nx*ny:nx*ny)
+	call inversetrans(fzs,nz,mxld,VerticalType)
+	x(j+(k-1)*nx:j+(k-1)*nx+(nz-1)*nx*ny:nx*ny)=fzs
+	enddo
+	enddo
+
+	do k=1,nz
+	do j=1,nx
+        fys=x(j+(k-1)*nx*ny:j+(ny-1)*nx+(k-1)*nx*ny:nx)
+	call inversetrans(fys,ny,mxl,HorizonType)
+        x(j+(k-1)*nx*ny:j+(ny-1)*nx+(k-1)*nx*ny:nx)=fys
+	enddo
+	enddo
+
+	do k=1,nz
+	do j=1,ny
+        fxs=x(1+(j-1)*nx+(k-1)*nx*ny:nx+(j-1)*nx+(k-1)*nx*ny)
+	call inversetrans(fxs,nx,mxl,HorizonType)
+        x(1+(j-1)*nx+(k-1)*nx*ny:nx+(j-1)*nx+(k-1)*nx*ny)=fxs
+	enddo
+	enddo
+	end subroutine

src/lsmrDataModule.f90

@@ -0,0 +1,24 @@
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! File lsmrDataModule.f90
+!
+! Defines real(dp) and a few constants for use in other modules.
+!
+! 24 Oct 2007: Allows floating-point precision dp to be defined
+!              in exactly one place (here).  Note that we need
+!                 use lsmrDataModule
+!              at the beginning of modules AND inside interfaces.
+!              zero and one are not currently used by LSMR,
+!              but this shows how they should be declared
+!              by a user routine that does need them.
+! 16 Jul 2010: LSMR version derived from LSQR equivalent.
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+module lsmrDataModule
+
+  implicit none
+
+  intrinsic                   ::      selected_real_kind
+  integer,  parameter, public :: dp = selected_real_kind(4)
+  real(dp), parameter, public :: zero = 0.0_dp, one = 1.0_dp
+
+end module lsmrDataModule

src/lsmrModule.f90

@@ -0,0 +1,754 @@
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! File lsmrModule.f90
+!
+!     LSMR
+!
+! LSMR   solves Ax = b or min ||Ax - b|| with or without damping,
+! using the iterative algorithm of David Fong and Michael Saunders:
+!     http://www.stanford.edu/group/SOL/software/lsmr.html
+!
+! Maintained by
+!     David Fong       <clfong@stanford.edu>
+!     Michael Saunders <saunders@stanford.edu>
+!     Systems Optimization Laboratory (SOL)
+!     Stanford University
+!     Stanford, CA 94305-4026, USA
+!
+! 17 Jul 2010: F90 LSMR derived from F90 LSQR and lsqr.m.
+! 07 Sep 2010: Local reorthogonalization now works (localSize > 0).
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+module lsmrModule
+
+  use  lsmrDataModule,    only : dp
+  use  lsmrblasInterface, only : dnrm2, dscal
+  implicit none
+  private
+  public   :: LSMR
+
+contains
+
+  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+!  subroutine LSMR  ( m, n, Aprod1, Aprod2, b, damp,                    &
+!                     atol, btol, conlim, itnlim, localSize, nout,      &
+!                     x, istop, itn, normA, condA, normr, normAr, normx )
+   subroutine LSMR  ( m, n, leniw, lenrw,iw,rw, b, damp,               &
+                     atol, btol, conlim, itnlim, localSize, nout,      &
+                     x, istop, itn, normA, condA, normr, normAr, normx )
+
+    integer,  intent(in) :: leniw
+    integer,  intent(in) :: lenrw
+    integer,  intent(in) :: iw(leniw)
+    real,     intent(in) :: rw(lenrw)
+
+    integer,  intent(in)  :: m, n, itnlim, localSize, nout
+    integer,  intent(out) :: istop, itn
+    real(dp), intent(in)  :: b(m)
+    real(dp), intent(out) :: x(n)
+    real(dp), intent(in)  :: atol, btol, conlim, damp
+    real(dp), intent(out) :: normA, condA, normr, normAr, normx
+
+    interface
+        subroutine aprod(mode,m,n,x,y,leniw,lenrw,iw,rw)          ! y := y + A*x
+         use lsmrDataModule, only : dp
+         integer,  intent(in)    :: mode,lenrw
+         integer,  intent(in)    :: leniw
+         real,     intent(in)    :: rw(lenrw)
+         integer,  intent(in)    :: iw(leniw)
+         integer,  intent(in)    :: m,n
+         real(dp), intent(inout)    :: x(n)
+         real(dp), intent(inout) :: y(m)
+       end subroutine aprod
+!      subroutine Aprod1(m,n,x,y)                   ! y := y + A*x
+!         use lsmrDataModule, only : dp
+!         integer,  intent(in)    :: m,n
+!         real(dp), intent(in)    :: x(n)
+!         real(dp), intent(inout) :: y(m)
+!       end subroutine Aprod1
+!
+!       subroutine Aprod2(m,n,x,y)                   ! x := x + A'*y
+!         use lsmrDataModule, only : dp
+!         integer,  intent(in)    :: m,n
+!         real(dp), intent(inout) :: x(n)
+!         real(dp), intent(in)    :: y(m)
+!       end subroutine Aprod2
+    end interface
+
+    !-------------------------------------------------------------------
+    ! LSMR  finds a solution x to the following problems:
+    !
+    ! 1. Unsymmetric equations:    Solve  A*x = b
+    !
+    ! 2. Linear least squares:     Solve  A*x = b
+    !                              in the least-squares sense
+    !
+    ! 3. Damped least squares:     Solve  (   A    )*x = ( b )
+    !                                     ( damp*I )     ( 0 )
+    !                              in the least-squares sense
+    !
+    ! where A is a matrix with m rows and n columns, b is an m-vector,
+    ! and damp is a scalar.  (All quantities are real.)
+    ! The matrix A is treated as a linear operator.  It is accessed
+    ! by means of subroutine calls with the following purpose:
+    !
+    ! call Aprod1(m,n,x,y)  must compute y = y + A*x  without altering x.
+    ! call Aprod2(m,n,x,y)  must compute x = x + A'*y without altering y.
+    !
+    ! LSMR uses an iterative method to approximate the solution.
+    ! The number of iterations required to reach a certain accuracy
+    ! depends strongly on the scaling of the problem.  Poor scaling of
+    ! the rows or columns of A should therefore be avoided where
+    ! possible.
+    !
+    ! For example, in problem 1 the solution is unaltered by
+    ! row-scaling.  If a row of A is very small or large compared to
+    ! the other rows of A, the corresponding row of ( A  b ) should be
+    ! scaled up or down.
+    !
+    ! In problems 1 and 2, the solution x is easily recovered
+    ! following column-scaling.  Unless better information is known,
+    ! the nonzero columns of A should be scaled so that they all have
+    ! the same Euclidean norm (e.g., 1.0).
+    !
+    ! In problem 3, there is no freedom to re-scale if damp is
+    ! nonzero.  However, the value of damp should be assigned only
+    ! after attention has been paid to the scaling of A.
+    !
+    ! The parameter damp is intended to help regularize
+    ! ill-conditioned systems, by preventing the true solution from
+    ! being very large.  Another aid to regularization is provided by
+    ! the parameter condA, which may be used to terminate iterations
+    ! before the computed solution becomes very large.
+    !
+    ! Note that x is not an input parameter.
+    ! If some initial estimate x0 is known and if damp = 0,
+    ! one could proceed as follows:
+    !
+    ! 1. Compute a residual vector     r0 = b - A*x0.
+    ! 2. Use LSMR to solve the system  A*dx = r0.
+    ! 3. Add the correction dx to obtain a final solution x = x0 + dx.
+    !
+    ! This requires that x0 be available before and after the call
+    ! to LSMR.  To judge the benefits, suppose LSMR takes k1 iterations
+    ! to solve A*x = b and k2 iterations to solve A*dx = r0.
+    ! If x0 is "good", norm(r0) will be smaller than norm(b).
+    ! If the same stopping tolerances atol and btol are used for each
+    ! system, k1 and k2 will be similar, but the final solution x0 + dx
+    ! should be more accurate.  The only way to reduce the total work
+    ! is to use a larger stopping tolerance for the second system.
+    ! If some value btol is suitable for A*x = b, the larger value
+    ! btol*norm(b)/norm(r0)  should be suitable for A*dx = r0.
+    !
+    ! Preconditioning is another way to reduce the number of iterations.
+    ! If it is possible to solve a related system M*x = b efficiently,
+    ! where M approximates A in some helpful way
+    ! (e.g. M - A has low rank or its elements are small relative to
+    ! those of A), LSMR may converge more rapidly on the system
+    !       A*M(inverse)*z = b,
+    ! after which x can be recovered by solving M*x = z.
+    !
+    ! NOTE: If A is symmetric, LSMR should not be used!
+    ! Alternatives are the symmetric conjugate-gradient method (CG)
+    ! and/or SYMMLQ.
+    ! SYMMLQ is an implementation of symmetric CG that applies to
+    ! any symmetric A and will converge more rapidly than LSMR.
+    ! If A is positive definite, there are other implementations of
+    ! symmetric CG that require slightly less work per iteration
+    ! than SYMMLQ (but will take the same number of iterations).
+    !
+    !
+    ! Notation
+    ! --------
+    ! The following quantities are used in discussing the subroutine
+    ! parameters:
+    !
+    ! Abar   =  (  A   ),        bbar  =  (b)
+    !           (damp*I)                  (0)
+    !
+    ! r      =  b - A*x,         rbar  =  bbar - Abar*x
+    !
+    ! normr  =  sqrt( norm(r)**2  +  damp**2 * norm(x)**2 )
+    !        =  norm( rbar )
+    !
+    ! eps    =  the relative precision of floating-point arithmetic.
+    !           On most machines, eps is about 1.0e-7 and 1.0e-16
+    !           in single and double precision respectively.
+    !           We expect eps to be about 1e-16 always.
+    !
+    ! LSMR  minimizes the function normr with respect to x.
+    !
+    !
+    ! Parameters
+    ! ----------
+    ! m       input      m, the number of rows in A.
+    !
+    ! n       input      n, the number of columns in A.
+    !
+    ! Aprod1, Aprod2     See above.
+    !
+    ! damp    input      The damping parameter for problem 3 above.
+    !                    (damp should be 0.0 for problems 1 and 2.)
+    !                    If the system A*x = b is incompatible, values
+    !                    of damp in the range 0 to sqrt(eps)*norm(A)
+    !                    will probably have a negligible effect.
+    !                    Larger values of damp will tend to decrease
+    !                    the norm of x and reduce the number of 
+    !                    iterations required by LSMR.
+    !
+    !                    The work per iteration and the storage needed
+    !                    by LSMR are the same for all values of damp.
+    !
+    ! b(m)    input      The rhs vector b.
+    !
+    ! x(n)    output     Returns the computed solution x.
+    !
+    ! atol    input      An estimate of the relative error in the data
+    !                    defining the matrix A.  For example, if A is
+    !                    accurate to about 6 digits, set atol = 1.0e-6.
+    !
+    ! btol    input      An estimate of the relative error in the data
+    !                    defining the rhs b.  For example, if b is
+    !                    accurate to about 6 digits, set btol = 1.0e-6.
+    !
+    ! conlim  input      An upper limit on cond(Abar), the apparent
+    !                    condition number of the matrix Abar.
+    !                    Iterations will be terminated if a computed
+    !                    estimate of cond(Abar) exceeds conlim.
+    !                    This is intended to prevent certain small or
+    !                    zero singular values of A or Abar from
+    !                    coming into effect and causing unwanted growth
+    !                    in the computed solution.
+    !
+    !                    conlim and damp may be used separately or
+    !                    together to regularize ill-conditioned systems.
+    !
+    !                    Normally, conlim should be in the range
+    !                    1000 to 1/eps.
+    !                    Suggested value:
+    !                    conlim = 1/(100*eps)  for compatible systems,
+    !                    conlim = 1/(10*sqrt(eps)) for least squares.
+    !
+    !         Note: Any or all of atol, btol, conlim may be set to zero.
+    !         The effect will be the same as the values eps, eps, 1/eps.
+    !
+    ! itnlim  input      An upper limit on the number of iterations.
+    !                    Suggested value:
+    !                    itnlim = n/2   for well-conditioned systems
+    !                                   with clustered singular values,
+    !                    itnlim = 4*n   otherwise.
+    !
+    ! localSize input    No. of vectors for local reorthogonalization.
+    !            0       No reorthogonalization is performed.
+    !           >0       This many n-vectors "v" (the most recent ones)
+    !                    are saved for reorthogonalizing the next v.
+    !                    localSize need not be more than min(m,n).
+    !                    At most min(m,n) vectors will be allocated.
+    !
+    ! nout    input      File number for printed output.  If positive,
+    !                    a summary will be printed on file nout.
+    !
+    ! istop   output     An integer giving the reason for termination:
+    !
+    !            0       x = 0  is the exact solution.
+    !                    No iterations were performed.
+    !
+    !            1       The equations A*x = b are probably compatible.
+    !                    Norm(A*x - b) is sufficiently small, given the
+    !                    values of atol and btol.
+    !
+    !            2       damp is zero.  The system A*x = b is probably
+    !                    not compatible.  A least-squares solution has
+    !                    been obtained that is sufficiently accurate,
+    !                    given the value of atol.
+    !
+    !            3       damp is nonzero.  A damped least-squares
+    !                    solution has been obtained that is sufficiently
+    !                    accurate, given the value of atol.
+    !
+    !            4       An estimate of cond(Abar) has exceeded conlim.
+    !                    The system A*x = b appears to be ill-conditioned,
+    !                    or there could be an error in Aprod1 or Aprod2.
+    !
+    !            5       The iteration limit itnlim was reached.
+    !
+    ! itn     output     The number of iterations performed.
+    !
+    ! normA   output     An estimate of the Frobenius norm of Abar.
+    !                    This is the square-root of the sum of squares
+    !                    of the elements of Abar.
+    !                    If damp is small and the columns of A
+    !                    have all been scaled to have length 1.0,
+    !                    normA should increase to roughly sqrt(n).
+    !                    A radically different value for normA may
+    !                    indicate an error in Aprod1 or Aprod2.
+    !
+    ! condA   output     An estimate of cond(Abar), the condition
+    !                    number of Abar.  A very high value of condA
+    !                    may again indicate an error in Aprod1 or Aprod2.
+    !
+    ! normr   output     An estimate of the final value of norm(rbar),
+    !                    the function being minimized (see notation
+    !                    above).  This will be small if A*x = b has
+    !                    a solution.
+    !
+    ! normAr  output     An estimate of the final value of
+    !                    norm( Abar'*rbar ), the norm of
+    !                    the residual for the normal equations.
+    !                    This should be small in all cases.  (normAr
+    !                    will often be smaller than the true value
+    !                    computed from the output vector x.)
+    !
+    ! normx   output     An estimate of norm(x) for the final solution x.
+    !
+    ! Subroutines and functions used              
+    ! ------------------------------
+    ! BLAS               dscal, dnrm2
+    ! USER               Aprod1, Aprod2
+    !
+    ! Precision
+    ! ---------
+    ! The number of iterations required by LSMR will decrease
+    ! if the computation is performed in higher precision.
+    ! At least 15-digit arithmetic should normally be used.
+    ! "real(dp)" declarations should normally be 8-byte words.
+    ! If this ever changes, the BLAS routines  dnrm2, dscal
+    ! (Lawson, et al., 1979) will also need to be changed.
+    !
+    !
+    ! Reference
+    ! ---------
+    ! http://www.stanford.edu/group/SOL/software/lsmr.html
+    ! ------------------------------------------------------------------
+    !
+    ! LSMR development:
+    ! 21 Sep 2007: Fortran 90 version of LSQR implemented.
+    !              Aprod1, Aprod2 implemented via f90 interface.
+    ! 17 Jul 2010: LSMR derived from LSQR and lsmr.m.
+    ! 07 Sep 2010: Local reorthogonalization now working.
+    !-------------------------------------------------------------------
+
+    intrinsic :: abs, dot_product, min, max, sqrt
+
+    ! Local arrays and variables
+    real(dp)  :: h(n), hbar(n), u(m), v(n), w(n), localV(n,min(localSize,m,n))
+    logical   :: damped, localOrtho, localVQueueFull, prnt, show
+    integer   :: i, localOrthoCount, localOrthoLimit, localPointer, localVecs, &
+                 pcount, pfreq
+    real(dp)  :: alpha, alphabar, alphahat, &
+                 beta, betaacute, betacheck, betad, betadd, betahat, &
+                 normb, c, cbar, chat, ctildeold, ctol,    &
+                 d, maxrbar, minrbar, normA2, &
+                 rho, rhobar, rhobarold, rhodold, rhoold, rhotemp, &
+                 rhotildeold, rtol, s, sbar, shat, stildeold, &
+                 t1, taud, tautildeold, test1, test2, test3, &
+                 thetabar, thetanew, thetatilde, thetatildeold, &
+                 zeta, zetabar, zetaold
+
+    ! Local constants
+    real(dp),         parameter :: zero = 0.0_dp, one = 1.0_dp
+    character(len=*), parameter :: enter = ' Enter LSMR.  '
+    character(len=*), parameter :: exitt = ' Exit  LSMR.  '
+    character(len=*), parameter :: msg(0:7) =                     &
+      (/ 'The exact solution is  x = 0                         ', &
+         'Ax - b is small enough, given atol, btol             ', &
+         'The least-squares solution is good enough, given atol', &
+         'The estimate of cond(Abar) has exceeded conlim       ', &
+         'Ax - b is small enough for this machine              ', &
+         'The LS solution is good enough for this machine      ', &
+         'Cond(Abar) seems to be too large for this machine    ', &
+         'The iteration limit has been reached                 ' /)
+    !-------------------------------------------------------------------
+
+
+    ! Initialize.
+
+    localVecs = min(localSize,m,n)
+    show      = nout > 0
+    if (show) then
+       write(nout, 1000) enter,m,n,damp,atol,conlim,btol,itnlim,localVecs
+    end if
+    
+    pfreq  = 20           ! print frequency (for repeating the heading)
+    pcount = 0            ! print counter
+    damped = damp > zero  !
+
+    !-------------------------------------------------------------------
+    ! Set up the first vectors u and v for the bidiagonalization.
+    ! These satisfy  beta*u = b,  alpha*v = A(transpose)*u.
+    !-------------------------------------------------------------------
+    u(1:m) = b(1:m)
+    v(1:n) = zero
+    x(1:n) = zero
+
+    alpha  = zero
+    beta   = dnrm2 (m, u, 1)
+
+    if (beta > zero) then
+       call dscal (m, (one/beta), u, 1)
+    !   call Aprod2(m, n, v, u)          ! v = A'*u
+        call aprod(2,m,n,v,u,leniw,lenrw,iw,rw)
+       alpha = dnrm2 (n, v, 1)
+    end if
+
+    if (alpha > zero) then
+       call dscal (n, (one/alpha), v, 1)
+       w = v
+    end if
+
+    normAr = alpha*beta
+    if (normAr == zero) go to 800
+
+    ! Initialization for local reorthogonalization.
+
+    localOrtho = .false.
+    if (localVecs > 0) then
+       localPointer    = 1
+       localOrtho      = .true.
+       localVQueueFull = .false.
+       localV(:,1)     = v
+    end if
+
+    ! Initialize variables for 1st iteration.
+
+    itn      = 0
+    zetabar  = alpha*beta
+    alphabar = alpha
+    rho      = 1
+    rhobar   = 1
+    cbar     = 1
+    sbar     = 0
+
+    h         = v
+    hbar(1:n) = zero
+    x(1:n)    = zero
+
+    ! Initialize variables for estimation of ||r||.
+
+    betadd      = beta
+    betad       = 0
+    rhodold     = 1
+    tautildeold = 0
+    thetatilde  = 0
+    zeta        = 0
+    d           = 0
+
+    ! Initialize variables for estimation of ||A|| and cond(A).
+
+    normA2  = alpha**2
+    maxrbar = 0_dp
+    minrbar = 1e+30_dp
+
+    ! Items for use in stopping rules.
+    normb  = beta
+    istop  = 0
+    ctol   = zero
+    if (conlim > zero) ctol = one/conlim
+    normr  = beta
+
+    ! Exit if b=0 or A'b = 0.
+
+    normAr = alpha * beta
+    if (normAr == 0) then
+       if (show) then
+          write(nout,'(a)') msg(1)
+       end if
+       return
+    end if
+
+    ! Heading for iteration log.
+
+    if (show) then
+       if (damped) then
+          write(nout,1300)
+       else
+          write(nout,1200)
+       end if
+       test1 = one
+       test2 = alpha/beta
+       write(nout, 1500) itn,x(1),normr,normAr,test1,test2
+    end if
+
+    !===================================================================
+    ! Main iteration loop.
+    !===================================================================
+    do
+       itn = itn + 1
+	  
+       !----------------------------------------------------------------
+       ! Perform the next step of the bidiagonalization to obtain the
+       ! next beta, u, alpha, v.  These satisfy
+       !     beta*u = A*v  - alpha*u,
+       !    alpha*v = A'*u -  beta*v.
+       !----------------------------------------------------------------
+       call dscal (m,(- alpha), u, 1)
+      ! call Aprod1(m, n, v, u)             ! u = A*v
+        call aprod ( 1,m,n,v,u,leniw,lenrw,iw,rw )      
+       beta   = dnrm2 (m, u, 1)
+
+       if (beta > zero) then
+          call dscal (m, (one/beta), u, 1)
+          if (localOrtho) then    ! Store v into the circular buffer localV.
+             call localVEnqueue   ! Store old v for local reorthog'n of new v.
+          end if
+          call dscal (n, (- beta), v, 1)
+
+          !call Aprod2(m, n, v, u)          ! v = A'*u
+          call aprod ( 2,m,n,v,u,leniw,lenrw,iw,rw )
+          if (localOrtho) then    ! Perform local reorthogonalization of V.
+             call localVOrtho     ! Local-reorthogonalization of new v.
+          end if
+          alpha  = dnrm2 (n, v, 1)
+          if (alpha > zero) then
+             call dscal (n, (one/alpha), v, 1)
+          end if
+       end if
+    
+       ! At this point, beta = beta_{k+1}, alpha = alpha_{k+1}.
+    
+       !----------------------------------------------------------------
+       ! Construct rotation Qhat_{k,2k+1}.
+
+       alphahat = d2norm(alphabar, damp)
+       chat     = alphabar/alphahat
+       shat     = damp/alphahat
+
+       ! Use a plane rotation (Q_i) to turn B_i to R_i.
+
+       rhoold   = rho
+       rho      = d2norm(alphahat, beta)
+       c        = alphahat/rho
+       s        = beta/rho
+       thetanew = s*alpha
+       alphabar = c*alpha
+
+       ! Use a plane rotation (Qbar_i) to turn R_i^T into R_i^bar.
+    
+       rhobarold = rhobar
+       zetaold   = zeta
+       thetabar  = sbar*rho
+       rhotemp   = cbar*rho
+       rhobar    = d2norm(cbar*rho, thetanew)
+       cbar      = cbar*rho/rhobar
+       sbar      = thetanew/rhobar
+       zeta      =   cbar*zetabar
+       zetabar   = - sbar*zetabar
+
+       ! Update h, h_hat, x.
+
+       hbar      = h - (thetabar*rho/(rhoold*rhobarold))*hbar
+       x         = x + (zeta/(rho*rhobar))*hbar
+       h         = v - (thetanew/rho)*h
+
+       ! Estimate ||r||.
+    
+       ! Apply rotation Qhat_{k,2k+1}.
+       betaacute =   chat* betadd
+       betacheck = - shat* betadd
+
+       ! Apply rotation Q_{k,k+1}.
+       betahat   =   c*betaacute
+       betadd    = - s*betaacute
+
+       ! Apply rotation Qtilde_{k-1}.
+       ! betad = betad_{k-1} here.
+
+       thetatildeold = thetatilde
+       rhotildeold   = d2norm(rhodold, thetabar)
+       ctildeold     = rhodold/rhotildeold
+       stildeold     = thetabar/rhotildeold
+       thetatilde    = stildeold* rhobar
+       rhodold       =   ctildeold* rhobar
+       betad         = - stildeold*betad + ctildeold*betahat
+
+       ! betad   = betad_k here.
+       ! rhodold = rhod_k  here.
+
+       tautildeold   = (zetaold - thetatildeold*tautildeold)/rhotildeold
+       taud          = (zeta - thetatilde*tautildeold)/rhodold
+       d             = d + betacheck**2
+       normr         = sqrt(d + (betad - taud)**2 + betadd**2)
+    
+       ! Estimate ||A||.
+       normA2        = normA2 + beta**2
+       normA         = sqrt(normA2)
+       normA2        = normA2 + alpha**2
+    
+       ! Estimate cond(A).
+       maxrbar       = max(maxrbar,rhobarold)
+       if (itn > 1) then 
+          minrbar    = min(minrbar,rhobarold)
+       end if
+       condA         = max(maxrbar,rhotemp)/min(minrbar,rhotemp)
+
+       !----------------------------------------------------------------
+       ! Test for convergence.
+       !----------------------------------------------------------------
+
+       ! Compute norms for convergence testing.
+       normAr  = abs(zetabar)
+       normx   = dnrm2(n, x, 1)
+
+       ! Now use these norms to estimate certain other quantities,
+       ! some of which will be small near a solution.
+
+       test1   = normr /normb
+       test2   = normAr/(normA*normr)
+       test3   =    one/condA
+       t1      =  test1/(one + normA*normx/normb)
+       rtol    =   btol + atol*normA*normx/normb
+
+       ! The following tests guard against extremely small values of
+       ! atol, btol or ctol.  (The user may have set any or all of
+       ! the parameters atol, btol, conlim  to 0.)
+       ! The effect is equivalent to the normAl tests using
+       ! atol = eps,  btol = eps,  conlim = 1/eps.
+
+       if (itn     >= itnlim) istop = 7
+       if (one+test3 <=  one) istop = 6
+       if (one+test2 <=  one) istop = 5
+       if (one+t1    <=  one) istop = 4
+
+       ! Allow for tolerances set by the user.
+
+       if (  test3   <= ctol) istop = 3
+       if (  test2   <= atol) istop = 2
+       if (  test1   <= rtol) istop = 1
+
+       !----------------------------------------------------------------
+       ! See if it is time to print something.
+       !----------------------------------------------------------------
+       prnt = .false.
+       if (show) then
+          if (n     <=        40) prnt = .true.
+          if (itn   <=        10) prnt = .true.
+          if (itn   >= itnlim-10) prnt = .true.
+          if (mod(itn,10)  ==  0) prnt = .true.
+          if (test3 <=  1.1*ctol) prnt = .true.
+          if (test2 <=  1.1*atol) prnt = .true.
+          if (test1 <=  1.1*rtol) prnt = .true.
+          if (istop /=         0) prnt = .true.
+
+          if (prnt) then        ! Print a line for this iteration
+             if (pcount >= pfreq) then  ! Print a heading first
+                pcount = 0
+                if (damped) then
+                   write(nout,1300)
+                else
+                   write(nout,1200)
+                end if
+             end if
+             pcount = pcount + 1
+             write(nout,1500) itn,x(1),normr,normAr,test1,test2,normA,condA
+          end if
+       end if
+
+       if (istop /= 0) exit
+    end do
+    !===================================================================
+    ! End of iteration loop.
+    !===================================================================
+
+    ! Come here if normAr = 0, or if normal exit.
+
+800 if (damped .and. istop==2) istop=3  ! Decide if istop = 2 or 3.
+    if (show) then                      ! Print the stopping condition.
+       write(nout, 2000)                &
+            exitt,istop,itn,            &
+            exitt,normA,condA,          &
+            exitt,normb, normx,         &
+            exitt,normr,normAr
+       write(nout, 3000)                &
+            exitt, msg(istop)
+    end if
+
+    return
+
+ 1000 format(// a, '     Least-squares solution of  Ax = b'       &
+      / ' The matrix  A  has', i7, ' rows   and', i7, ' columns'  &
+      / ' damp   =', es22.14                                      &
+      / ' atol   =', es10.2, 15x,        'conlim =', es10.2       &
+      / ' btol   =', es10.2, 15x,        'itnlim =', i10          &
+      / ' localSize (no. of vectors for local reorthogonalization) =', i7)
+ 1200 format(/ "   Itn       x(1)            norm r         A'r   ", &
+      ' Compatible    LS      norm A    cond A')
+ 1300 format(/ "   Itn       x(1)           norm rbar    Abar'rbar", &
+      ' Compatible    LS    norm Abar cond Abar')
+ 1500 format(i6, 2es17.9, 5es10.2)
+ 2000 format(/ a, 5x, 'istop  =', i2,   15x, 'itn    =', i8      &
+      /      a, 5x, 'normA  =', es12.5, 5x, 'condA  =', es12.5   &
+      /      a, 5x, 'normb  =', es12.5, 5x, 'normx  =', es12.5   &
+      /      a, 5x, 'normr  =', es12.5, 5x, 'normAr =', es12.5)
+ 3000 format(a, 5x, a)
+
+  contains
+
+    function d2norm( a, b )
+
+      real(dp)             :: d2norm
+      real(dp), intent(in) :: a, b
+
+      !-------------------------------------------------------------------
+      ! d2norm returns sqrt( a**2 + b**2 )
+      ! with precautions to avoid overflow.
+      !
+      ! 21 Mar 1990: First version.
+      ! 17 Sep 2007: Fortran 90 version.
+      ! 24 Oct 2007: User real(dp) instead of compiler option -r8.
+      !-------------------------------------------------------------------
+
+      intrinsic            :: abs, sqrt
+      real(dp)             :: scale
+      real(dp), parameter  :: zero = 0.0_dp
+
+      scale = abs(a) + abs(b)
+      if (scale == zero) then
+         d2norm = zero
+      else
+         d2norm = scale*sqrt((a/scale)**2 + (b/scale)**2)
+      end if
+
+    end function d2norm
+
+    !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    subroutine localVEnqueue
+
+      ! Store v into the circular buffer localV.
+
+      if (localPointer < localVecs) then
+         localPointer = localPointer + 1
+      else
+         localPointer = 1
+         localVQueueFull = .true.
+      end if
+      localV(:,localPointer) = v
+
+    end subroutine localVEnqueue
+
+    !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+    subroutine localVOrtho
+
+      ! Perform local reorthogonalization of current v.
+
+      real(dp)  :: d
+
+      if (localVQueueFull) then
+         localOrthoLimit = localVecs
+      else
+         localOrthoLimit = localPointer
+      end if
+
+      do localOrthoCount = 1, localOrthoLimit
+         d = dot_product(v,localV(:,localOrthoCount))
+         v = v    -    d * localV(:,localOrthoCount)
+      end do
+
+    end subroutine localVOrtho
+
+  end subroutine LSMR
+
+  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+end module LSMRmodule

src/lsmrblas.f90

@@ -0,0 +1,360 @@
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!     File lsmrblas.f90   (double precision)
+!
+!     This file contains the following BLAS routines
+!        dcopy, ddot, dnrm2, dscal
+!     required by subroutines LSMR and Acheck.
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!
+!! DCOPY copies a vector X to a vector Y.
+!
+!  Discussion:
+!    This routine uses double precision real arithmetic.
+!    The routine uses unrolled loops for increments equal to one.
+!
+!  Modified:
+!    16 May 2005
+!
+!  Author:
+!    Jack Dongarra
+!    Fortran90 translation by John Burkardt.
+!
+!  Reference:
+!
+!    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
+!    LINPACK User's Guide,
+!    SIAM, 1979,
+!    ISBN13: 978-0-898711-72-1,
+!    LC: QA214.L56.
+!
+!    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
+!    Algorithm 539, 
+!    Basic Linear Algebra Subprograms for Fortran Usage,
+!    ACM Transactions on Mathematical Software, 
+!    Volume 5, Number 3, September 1979, pages 308-323.
+!
+!  Parameters:
+!
+!    Input, integer N, the number of elements in DX and DY.
+!
+!    Input, real ( kind = 8 ) DX(*), the first vector.
+!
+!    Input, integer INCX, the increment between successive entries of DX.
+!
+!    Output, real ( kind = 8 ) DY(*), the second vector.
+!
+!    Input, integer INCY, the increment between successive entries of DY.
+
+
+      subroutine  dcopy(n,dx,incx,dy,incy)
+
+      implicit none
+!      double precision dx(*),dy(*)
+      real(4) dx(*),dy(*)
+      integer i,incx,incy,ix,iy,m,n
+
+      if ( n <= 0 ) then
+         return
+      end if
+
+      if ( incx == 1 .and. incy == 1 ) then
+
+         m = mod ( n, 7 )
+
+         if ( m /= 0 ) then
+            dy(1:m) = dx(1:m)
+         end if
+
+         do i = m+1, n, 7
+            dy(i) = dx(i)
+            dy(i + 1) = dx(i + 1)
+            dy(i + 2) = dx(i + 2)
+            dy(i + 3) = dx(i + 3)
+            dy(i + 4) = dx(i + 4)
+            dy(i + 5) = dx(i + 5)
+            dy(i + 6) = dx(i + 6)
+         end do
+
+        else
+
+           if ( 0 <= incx ) then
+              ix = 1
+           else
+              ix = ( -n + 1 ) * incx + 1
+           end if
+
+           if ( 0 <= incy ) then
+              iy = 1
+           else
+              iy = ( -n + 1 ) * incy + 1
+           end if
+
+           do i = 1, n
+              dy(iy) = dx(ix)
+              ix = ix + incx
+              iy = iy + incy
+           end do
+        end if
+        return
+end subroutine dcopy
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!! DDOT forms the dot product of two vectors.
+!
+!  Discussion:
+!    This routine uses double precision real arithmetic.
+!    This routine uses unrolled loops for increments equal to one.
+!
+!  Modified:
+!    16 May 2005
+!
+!  Author:
+!    Jack Dongarra
+!    Fortran90 translation by John Burkardt.
+!
+!  Reference:
+!    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
+!    LINPACK User's Guide,
+!    SIAM, 1979,
+!    ISBN13: 978-0-898711-72-1,
+!    LC: QA214.L56.
+!
+!    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
+!    Algorithm 539, 
+!    Basic Linear Algebra Subprograms for Fortran Usage,
+!    ACM Transactions on Mathematical Software, 
+!    Volume 5, Number 3, September 1979, pages 308-323.
+!
+!  Parameters:
+!
+!    Input, integer N, the number of entries in the vectors.
+!
+!    Input, real ( kind = 8 ) DX(*), the first vector.
+!
+!    Input, integer INCX, the increment between successive entries in DX.
+!
+!    Input, real ( kind = 8 ) DY(*), the second vector.
+!
+!    Input, integer INCY, the increment between successive entries in DY.
+!
+!    Output, real ( kind = 8 ) DDOT, the sum of the product of the 
+!    corresponding entries of DX and DY.
+
+
+     ! double precision function ddot(n,dx,incx,dy,incy)
+      real(4) function ddot(n,dx,incx,dy,incy)
+
+      implicit         none
+    ! double precision dx(*),dy(*),dtemp
+      real(4) dx(*),dy(*),dtemp
+      integer          i,incx,incy,ix,iy,m,n
+
+      ddot = 0.0d0
+      dtemp = 0.0d0
+      if ( n <= 0 ) then
+         return
+      end if
+
+!  Code for unequal increments or equal increments
+!  not equal to 1.
+
+      if ( incx /= 1 .or. incy /= 1 ) then
+
+         if ( 0 <= incx ) then
+            ix = 1
+         else
+            ix = ( - n + 1 ) * incx + 1
+         end if
+
+         if ( 0 <= incy ) then
+            iy = 1
+         else
+            iy = ( - n + 1 ) * incy + 1
+         end if
+
+         do i = 1, n
+            dtemp = dtemp + dx(ix) * dy(iy)
+            ix = ix + incx
+            iy = iy + incy
+         end do
+
+!  Code for both increments equal to 1.
+
+        else
+
+           m = mod ( n, 5 )
+
+           do i = 1, m
+              dtemp = dtemp + dx(i) * dy(i)
+           end do
+
+           do i = m+1, n, 5
+              dtemp = dtemp + dx(i)*dy(i) + dx(i+1)*dy(i+1) + dx(i+2)*dy(i+2) &
+                                          + dx(i+3)*dy(i+3) + dx(i+4)*dy(i+4)
+           end do
+
+        end if
+
+        ddot = dtemp
+        return
+end function ddot
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!*****************************************************************************80
+!
+!! DNRM2 returns the euclidean norm of a vector.
+!
+!  Discussion:
+!    This routine uses double precision real arithmetic.
+!     DNRM2 ( X ) = sqrt ( X' * X )
+!
+!  Modified:
+!    16 May 2005
+!
+!  Author:
+!    Sven Hammarling
+!    Fortran90 translation by John Burkardt.
+!
+!  Reference:
+!    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
+!    LINPACK User's Guide,
+!    SIAM, 1979,
+!    ISBN13: 978-0-898711-72-1,
+!    LC: QA214.L56.
+!
+!    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
+!    Algorithm 539, 
+!    Basic Linear Algebra Subprograms for Fortran Usage,
+!    ACM Transactions on Mathematical Software,
+!    Volume 5, Number 3, September 1979, pages 308-323.
+!
+!  Parameters:
+!
+!    Input, integer N, the number of entries in the vector.
+!
+!    Input, real ( kind = 8 ) X(*), the vector whose norm is to be computed.
+!
+!    Input, integer INCX, the increment between successive entries of X.
+!
+!    Output, real ( kind = 8 ) DNRM2, the Euclidean norm of X.
+!
+
+    !  double precision function dnrm2 ( n, x, incx)
+      real(4) function dnrm2 ( n, x, incx)
+      implicit         none
+      integer          ix,n,incx
+     ! double precision x(*), ssq,absxi,norm,scale
+      real(4) x(*), ssq,absxi,norm,scale
+
+      if ( n < 1 .or. incx < 1 ) then
+         norm  = 0.d0
+      else if ( n == 1 ) then
+         norm  = abs ( x(1) )
+      else
+         scale = 0.d0
+         ssq = 1.d0
+
+         do ix = 1, 1 + ( n - 1 )*incx, incx
+            if ( x(ix) /= 0.d0 ) then
+               absxi = abs ( x(ix) )
+               if ( scale < absxi ) then
+                  ssq = 1.d0 + ssq * ( scale / absxi )**2
+                  scale = absxi
+               else
+                  ssq = ssq + ( absxi / scale )**2
+               end if
+            end if
+         end do
+         norm  = scale * sqrt ( ssq )
+      end if
+
+      dnrm2 = norm
+      return
+end function dnrm2
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! DSCAL scales a vector by a constant.
+!
+!  Discussion:
+!    This routine uses double precision real arithmetic.
+!
+!  Modified:
+!    08 April 1999
+!
+!  Author:
+!    Jack Dongarra
+!    Fortran90 translation by John Burkardt.
+!
+!  Reference:
+!    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
+!    LINPACK User's Guide,
+!    SIAM, 1979,
+!    ISBN13: 978-0-898711-72-1,
+!    LC: QA214.L56.
+!
+!    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
+!    Algorithm 539, 
+!    Basic Linear Algebra Subprograms for Fortran Usage,
+!    ACM Transactions on Mathematical Software,
+!    Volume 5, Number 3, September 1979, pages 308-323.
+!
+!  Parameters:
+!
+!    Input, integer N, the number of entries in the vector.
+!
+!    Input, real ( kind = 8 ) SA, the multiplier.
+!
+!    Input/output, real ( kind = 8 ) X(*), the vector to be scaled.
+!
+!    Input, integer INCX, the increment between successive entries of X.
+!
+
+      subroutine  dscal(n,sa,x,incx)
+
+      implicit none
+
+      integer i
+      integer incx
+      integer ix
+      integer m
+      integer n
+      !double precision sa
+      !double precision x(*)
+
+      real(4) sa
+      real(4) x(*)
+      if ( n <= 0 ) then
+         return
+      else if ( incx == 1 ) then
+         m = mod ( n, 5 )
+         x(1:m) = sa * x(1:m)
+
+         do i = m+1, n, 5
+            x(i)   = sa * x(i)
+            x(i+1) = sa * x(i+1)
+            x(i+2) = sa * x(i+2)
+            x(i+3) = sa * x(i+3)
+            x(i+4) = sa * x(i+4)
+         end do
+      else
+         if ( 0 <= incx ) then
+            ix = 1
+         else
+            ix = ( - n + 1 ) * incx + 1
+         end if
+
+         do i = 1, n
+            x(ix) = sa * x(ix)
+            ix = ix + incx
+         end do
+
+      end if
+
+      return
+end subroutine dscal
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

src/lsmrblasInterface.f90

@@ -0,0 +1,41 @@
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! File lsmrblasInterface.f90
+!
+!    BLAS1 Interfaces:   ddot    dnrm2    dscal
+!
+! Maintained by Michael Saunders <saunders@stanford.edu>.
+!
+! 19 Dec 2008: lsqrblasInterface module implemented.
+!              Metcalf and Reid recommend putting interfaces in a module.
+! 16 Jul 2010: LSMR version derived from LSQR equivalent.
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+module lsmrblasInterface
+
+  implicit none
+  public   :: ddot, dnrm2, dscal
+
+  interface                              ! Level 1 BLAS
+     function ddot  (n,dx,incx,dy,incy)
+       use lsmrDataModule, only : dp
+       integer,  intent(in)    :: n,incx,incy
+       real(dp), intent(in)    :: dx(*),dy(*)
+       real(dp)                :: ddot
+     end function ddot
+
+     function dnrm2 (n,dx,incx)
+       use lsmrDataModule, only : dp
+       integer,  intent(in)    :: n,incx
+       real(dp), intent(in)    :: dx(*)
+       real(dp)                :: dnrm2
+     end function dnrm2
+
+     subroutine dscal (n,sa,x,incx)
+       use lsmrDataModule, only : dp
+       integer,  intent(in)    :: n,incx
+       real(dp), intent(in)    :: sa
+       real(dp), intent(inout) :: x(*)
+     end subroutine dscal
+  end interface
+
+end module lsmrblasInterface

src/main.f90

@@ -0,0 +1,756 @@
+      ! CODE FOR SURFACE WAVE TOMOGRAPHY USING DISPERSION MEASUREMENTS
+      ! VERSION:
+      !      1.0 
+      ! AUTHOR:
+      !      HONGJIAN FANG. fanghj@mail.ustc.edu.cn
+      ! PURPOSE:
+      !      DIRECTLY INVERT SURFACE WAVE DISPERSION MEASUREMENTS FOR 3-D
+      ! STUCTURE WITHOUT THE INTERMEDIATE STEP OF CONSTUCTION THE PHASE
+      ! OR GROUP VELOCITY MAPS.
+      ! REFERENCE:
+      ! 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.
+      ! HISTORY:
+      !       2015/01/31 START TO REORGONIZE THE MESSY CODE   
+      !
+
+        program SurfTomo
+        use lsmrModule, only:lsmr
+	use  lsmrblasInterface, only : dnrm2
+        implicit none
+
+! VARIABLE DEFINE
+
+            character inputfile*80
+            character logfile*100
+            character outmodel*100
+            character outsyn*100
+            logical ex
+            character dummy*40
+            character datafile*80
+
+            integer nx,ny,nz
+            real goxd,gozd
+            real dvxd,dvzd
+            integer nsrc,nrc
+            real weight,weight0
+            real damp
+            real minthk
+            integer kmax,kmaxRc,kmaxRg,kmaxLc,kmaxLg
+            real*8,dimension(:),allocatable:: tRc,tRg,tLc,tLg
+            real,dimension(:),allocatable:: depz
+            integer itn
+            integer nout
+            integer localSize
+            real mean,std_devs,balances,balanceb
+            integer msurf
+            integer maxlevel,maxleveld
+            real,parameter:: tolr=1e-4
+            real,dimension(:),allocatable:: obst,dsyn,cbst,wt,dtres,dist,datweight
+         	real,dimension(:),allocatable:: pvall,depRp,pvRp
+        	real sta1_lat,sta1_lon,sta2_lat,sta2_lon
+        	real dist1
+                integer dall
+         	integer istep
+         	real,parameter :: pi=3.1415926535898
+         	integer checkstat
+         	integer ii,jj,kk
+                real, dimension (:,:), allocatable :: scxf,sczf
+                real, dimension (:,:,:), allocatable :: rcxf,rczf
+	        integer,dimension(:,:),allocatable::wavetype,igrt,nrc1
+         	integer,dimension(:),allocatable::nsrc1,knum1
+         	integer,dimension(:,:),allocatable::periods
+         	real,dimension(:),allocatable::rw
+         	integer,dimension(:),allocatable::iw,col
+                real,dimension(:),allocatable::dv,norm
+                real,dimension(:,:,:),allocatable::vsf
+                real,dimension(:,:,:),allocatable::vsftrue
+        	character strf
+        	integer veltp,wavetp
+        	real velvalue
+        	integer knum,knumo,err
+        	integer istep1,istep2
+        	integer period
+        	integer knumi,srcnum,count1
+        	integer HorizonType,VerticalType
+        	character line*200
+            integer iter,maxiter
+            integer iiter,initer
+            integer maxnar
+            real acond
+            real anorm
+            real arnorm
+            real rnorm
+            real xnorm
+            character str1
+            real atol,btol
+            real conlim
+            integer istop
+            integer itnlim
+            integer lenrw,leniw
+            integer nar,nar_tmp,nars
+            integer count3,nvz,nvx
+            integer m,maxvp,n
+            integer i,j,k
+            real Minvel,MaxVel
+            real spfra
+            integer domain,normt
+            real noiselevel
+            integer ifsyn
+            integer writepath
+            real averdws
+            real maxnorm
+	    real threshold,threshold0
+
+! OPEN FILES FIRST TO OUTPUT THE PROCESS
+            open(34,file='IterVel.out')
+            nout=36
+            open(nout,file='lsmr.txt')
+
+! OUTPUT PROGRAM INFOMATION            
+             write(*,*)
+             write(*,*),'                         S U R F  T O M O'
+             write(*,*),'PLEASE contact Hongjain Fang &
+                 (fanghj@mail.ustc.edu.cn) if you find any bug'
+             write(*,*)
+
+! READ INPUT FILE
+            if (iargc() < 1) then
+                write(*,*) 'input file [SurfTomo.in(default)]:'
+                read(*,'(a)') inputfile
+                if (len_trim(inputfile) <=1 ) then
+                    inputfile = 'SurfTomo.in'
+                else
+                    inputfile = inputfile(1:len_trim(inputfile))
+                endif
+            else
+                call getarg(1,inputfile)
+            endif
+            inquire(file = inputfile, exist = ex)
+            if (.not. ex) stop 'unable to open the inputfile'
+
+            open(10,file=inputfile,status='old')
+            read(10,'(a30)')dummy
+            read(10,'(a30)')dummy
+            read(10,'(a30)')dummy
+            read(10,*)datafile
+            read(10,*) nx,ny,nz
+            read(10,*) goxd,gozd
+            read(10,*) dvxd,dvzd
+            read(10,*) nsrc
+            read(10,*) weight0,damp
+            read(10,*) minthk
+            read(10,*) Minvel,Maxvel
+            read(10,*) domain,normt
+	        read(10,*) HorizonType,VerticalType
+            read(10,*) maxlevel,maxleveld
+            read(10,*) maxiter,initer
+            read(10,*) spfra
+            read(10,*) kmaxRc
+            write(*,*)  'model origin:latitude,longitue'
+            write(*,'(2f10.4)') goxd,gozd
+            write(*,*) 'grid spacing:latitude,longitue'
+            write(*,'(2f10.4)') dvxd,dvzd
+            write(*,*) 'model dimension:nx,ny,nz'
+            write(*,'(3i5)') nx,ny,nz
+	        write(logfile,'(a,a)')trim(inputfile),'.log' 
+            open(66,file=logfile)
+            write(66,*)
+            write(66,*),'                         S U R F  T O M O'
+            write(66,*),'PLEASE contact Hongjain Fang &
+                (fanghj@mail.ustc.edu.cn) if you find any bug'
+            write(66,*)
+            write(66,*) 'model origin:latitude,longitue'
+            write(66,'(2f10.4)') goxd,gozd
+            write(66,*) 'grid spacing:latitude,longitue'
+            write(66,'(2f10.4)') dvxd,dvzd
+            write(66,*) 'model dimension:nx,ny,nz'
+            write(66,'(3i5)') nx,ny,nz
+            if(kmaxRc.gt.0)then
+               allocate(tRc(kmaxRc),&
+                stat=checkstat)
+               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(66,*)'Rayleigh wave phase velocity used,periods:(s)'
+            write(66,'(50f6.2)')(tRc(i),i=1,kmaxRc)
+            endif
+            read(10,*)kmaxRg
+            if(kmaxRg.gt.0)then
+               allocate(tRg(kmaxRg), stat=checkstat)
+               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(66,*)'Rayleigh wave group velocity used,periods:(s)'
+            write(66,'(50f6.2)')(tRg(i),i=1,kmaxRg)
+            endif
+            read(10,*)kmaxLc
+            if(kmaxLc.gt.0)then
+               allocate(tLc(kmaxLc), stat=checkstat)
+               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(66,*)'Love wave phase velocity used,periods:(s)'
+            write(66,'(50f6.2)')(tLc(i),i=1,kmaxLc)
+            endif
+            read(10,*)kmaxLg
+            if(kmaxLg.gt.0)then
+               allocate(tLg(kmaxLg), stat=checkstat)
+               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(66,*)'Love wave group velocity used,periods:(s)'
+            write(66,'(50f6.2)')(tLg(i),i=1,kmaxLg)
+            endif
+            read(10,*)ifsyn
+            read(10,*)noiselevel
+	    read(10,*) threshold0
+            close(10)
+            nrc=nsrc
+            kmax=kmaxRc+kmaxRg+kmaxLc+kmaxLg
+
+! READ MEASUREMENTS            
+            open(unit=87,file=datafile,status='old')
+            allocate(scxf(nsrc,kmax),sczf(nsrc,kmax),&
+            rcxf(nrc,nsrc,kmax),rczf(nrc,nsrc,kmax),stat=checkstat)
+            if(checkstat > 0)then
+            write(6,*)'error with allocate'
+            endif
+            allocate(periods(nsrc,kmax),wavetype(nsrc,kmax),&
+            nrc1(nsrc,kmax),nsrc1(kmax),knum1(kmax),&
+            igrt(nsrc,kmax),stat=checkstat)
+            if(checkstat > 0)then
+            write(6,*)'error with allocate'
+            endif
+            allocate(obst(nrc*nsrc*kmax),dist(nrc*nsrc*kmax),&
+            stat=checkstat)
+            allocate(pvall(nrc*nsrc*kmax),depRp(nrc*nsrc*kmax),&
+            pvRp(nrc*nsrc*kmax),stat=checkstat)
+            IF(checkstat > 0)THEN
+            write(6,*)'error with allocate'
+            ENDIF
+            istep=0
+            istep2=0
+            dall=0
+            knumo=12345
+	        knum=0
+	        istep1=0
+            do 
+            read(87,'(a)',iostat=err) line
+            if(err.eq.0) then
+            if(line(1:1).eq.'#') then
+            read(line,*) str1,sta1_lat,sta1_lon,period,wavetp,veltp
+            if(wavetp.eq.2.and.veltp.eq.0) knum=period
+            if(wavetp.eq.2.and.veltp.eq.1) knum=kmaxRc+period
+            if(wavetp.eq.1.and.veltp.eq.0) knum=kmaxRg+kmaxRc+period
+            if(wavetp.eq.1.and.veltp.eq.1) knum=kmaxLc+kmaxRg+&
+               kmaxRc+period
+            if(knum.ne.knumo) then
+            istep=0
+            istep2=istep2+1
+            endif
+            istep=istep+1
+            istep1=0
+            sta1_lat=(90.0-sta1_lat)*pi/180.0
+            sta1_lon=sta1_lon*pi/180.0
+            scxf(istep,knum)=sta1_lat
+            sczf(istep,knum)=sta1_lon
+            periods(istep,knum)=period
+            wavetype(istep,knum)=wavetp
+            igrt(istep,knum)=veltp
+            nsrc1(knum)=istep
+            knum1(istep2)=knum
+            knumo=knum
+            else
+            read(line,*) sta2_lat,sta2_lon,velvalue
+            istep1=istep1+1
+            dall=dall+1
+            sta2_lat=(90.0-sta2_lat)*pi/180.0
+            sta2_lon=sta2_lon*pi/180.0
+            rcxf(istep1,istep,knum)=sta2_lat
+            rczf(istep1,istep,knum)=sta2_lon
+            call delsph(sta1_lat,sta1_lon,sta2_lat,sta2_lon,dist1)
+	    dist(dall)=dist1
+            obst(dall)=dist1/velvalue
+            pvall(dall)=velvalue
+            nrc1(istep,knum)=istep1
+            endif
+            else
+            exit
+            endif
+            enddo
+            close(87)
+            allocate(depz(nz), stat=checkstat)
+            maxnar = dall*nx*ny*nz*spfra!sparsity fraction
+            maxvp = (nx-2)*(ny-2)*(nz-1)
+            allocate(dv(maxvp), stat=checkstat)
+            allocate(norm(maxvp), stat=checkstat)
+            allocate(vsf(nx,ny,nz), stat=checkstat)
+            allocate(vsftrue(nx,ny,nz), stat=checkstat)
+        	
+        	allocate(rw(maxnar), stat=checkstat)
+        	if(checkstat > 0)then
+        	   write(6,*)'error with allocate: real rw'
+        	endif
+        	allocate(iw(2*maxnar+1), stat=checkstat)
+        	if(checkstat > 0)then
+        	   write(6,*)'error with allocate: integer iw'
+        	endif
+        	allocate(col(maxnar), stat=checkstat)
+        	if(checkstat > 0)then
+        	   write(6,*)'error with allocate:  integer iw'
+        	endif
+           allocate(cbst(dall+maxvp),dsyn(dall),datweight(dall),wt(dall+maxvp),dtres(dall+maxvp),&
+                stat=checkstat)
+
+! MEASUREMENTS STATISTICS AND READ INITIAL MODEL           
+            write(*,'(a,i7)') 'Number of all measurements',dall
+
+            open(10,file='MOD',status='old')
+            read(10,*) (depz(i),i=1,nz)
+            do k = 1,nz
+            do j = 1,ny
+            read(10,*)(vsf(i,j,k),i=1,nx)
+            enddo
+            enddo
+            close(10)
+            write(*,*) 'grid points in depth direction:(km)'
+            write(*,'(50f6.2)') depz
+
+
+
+! CHECKERBOARD TEST
+            if (ifsyn == 1) then
+            write(*,*) 'Checkerboard Resolution Test Begin'
+            vsftrue = vsf
+
+            open(11,file='MOD.true',status='old')
+            do k = 1,nz
+            do j = 1,ny
+            read(11,*) (vsftrue(i,j,k),i=1,nx)
+            enddo
+            enddo
+            close(11)
+
+            call synthetic(nx,ny,nz,maxvp,vsftrue,obst,&
+                goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
+                tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
+                scxf,sczf,rcxf,rczf,nrc1,nsrc1,knum1,kmax,&
+                nsrc,nrc,noiselevel)
+            endif
+
+
+
+! ITERATE UNTILL CONVERGE
+            writepath = 0
+            do iter = 1,maxiter
+            iw = 0
+            rw = 0.0
+            col = 0
+
+! COMPUTE SENSITIVITY MATRIX
+            if (iter == maxiter) then
+                writepath = 1
+                open(40,file='raypath.out')
+            endif
+            write(*,*) 'computing sensitivity matrix...'
+            call CalSurfG(nx,ny,nz,maxvp,vsf,iw,rw,col,dsyn,&
+                goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
+                tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
+                scxf,sczf,rcxf,rczf,nrc1,nsrc1,knum1,kmax,&
+                nsrc,nrc,nar,domain,&
+                maxlevel,maxleveld,HorizonType,VerticalType,writepath)
+
+	        do i = 1,dall
+	        cbst(i) = obst(i) - dsyn(i)
+	        enddo
+
+		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
+		cbst(i) = cbst(i)*datweight(i)
+		enddo
+
+		do i = 1,nar
+		rw(i) = rw(i)*datweight(iw(1+i))
+		enddo
+ 
+            	norm=0
+            	do i=1,nar
+            	norm(col(i))=norm(col(i))+abs(rw(i))
+            	enddo
+            	averdws=0
+            	maxnorm=0
+            	do i=1,maxvp
+            	averdws =  averdws+norm(i)
+            	if(norm(i)>maxnorm) maxnorm=norm(i)
+            	enddo
+            	averdws=averdws/maxvp
+            	write(66,*)'Maximum and Average DWS values:',maxnorm,averdws
+            	write(66,*)'Threshold is:',threshold
+            
+! WRITE OUT RESIDUAL FOR THE FIRST AND LAST ITERATION
+               if(iter.eq.1) then
+               open(88,file='residualFirst.dat')
+               do i=1,dall
+               write(88,*) dist(i),dsyn(i),obst(i), &
+               dsyn(i)*datweight(i),obst(i)*datweight(i),datweight(i)
+               enddo
+               close(88)
+               endif
+               if(iter.eq.maxiter) then
+               open(88,file='residualLast.dat')
+               do i=1,dall
+               write(88,*) dist(i),dsyn(i),obst(i), &
+               dsyn(i)*datweight(i),obst(i)*datweight(i),datweight(i)
+               enddo
+               close(88)
+               endif
+ 
+         
+! ADDING REGULARIZATION TERM
+       	       weight=dnrm2(dall,cbst,1)**2/dall*weight0
+               nar_tmp=nar
+     	       nars=0
+             !  if (domain == 0 .and. normt==0) then
+               if (domain == 0) then
+               do i=1,maxvp
+               rw(nar+i)=weight
+               iw(1+nar+i)=dall+i
+               col(nar+i)=i
+               cbst(dall+i)=0
+               enddo
+               nar = nar + maxvp
+                m = dall + maxvp
+                n = maxvp
+             !	elseif(domain == 0 .and. normt/=0) then
+             !	do i=1,maxvp
+             !   rw(nar+i)=weight
+             !   iw(1+nar+i)=dall+i
+             !   col(nar+i)=i
+             !   cbst(dall+i)=0
+             !   enddo
+             !   nar = nar + maxvp
+             !    m = dall + maxvp
+             !    n = maxvp
+                 else
+                 count3=0
+                 nvz=ny-2
+                 nvx=nx-2
+                 do k=1,nz-1
+                 do j=1,nvz
+                 do i=1,nvx
+                 if(i==1.or.i==nvx.or.j==1.or.j==nvz.or.k==1.or.k==nz-1)then
+                 count3=count3+1
+                 col(nar+1)=(k-1)*nvz*nvx+(j-1)*nvx+i
+                 rw(nar+1)=2.0*weight
+                 iw(1+nar+1)=dall+count3
+                 cbst(dall+count3)=0
+                 nar=nar+1
+                 else
+                 count3=count3+1
+                 col(nar+1)=(k-1)*nvz*nvx+(j-1)*nvx+i
+                 rw(nar+1)=6.0*weight
+                 iw(1+nar+1)=dall+count3
+                 rw(nar+2)=-1.0*weight
+                 iw(1+nar+2)=dall+count3
+                 col(nar+2)=(k-1)*nvz*nvx+(j-1)*nvx+i-1
+                 rw(nar+3)=-1.0*weight
+                 iw(1+nar+3)=dall+count3
+                 col(nar+3)=(k-1)*nvz*nvx+(j-1)*nvx+i+1
+                 rw(nar+4)=-1.0*weight
+                 iw(1+nar+4)=dall+count3
+                 col(nar+4)=(k-1)*nvz*nvx+(j-2)*nvx+i
+                 rw(nar+5)=-1.0*weight
+                 iw(1+nar+5)=dall+count3
+                 col(nar+5)=(k-1)*nvz*nvx+j*nvx+i
+                 rw(nar+6)=-1.0*weight
+                 iw(1+nar+6)=dall+count3
+                 col(nar+6)=(k-2)*nvz*nvx+(j-1)*nvx+i
+                 rw(nar+7)=-1.0*weight
+                 iw(1+nar+7)=dall+count3
+                 col(nar+7)=k*nvz*nvx+(j-1)*nvx+i
+                 cbst(dall+count3)=0
+                 nar=nar+7
+                 endif
+                 enddo
+                 enddo
+                 enddo
+                   m = dall + count3
+                   n = maxvp
+                   nars = nar - nar_tmp
+                   rw(nar+1:nar+nars) = rw(nar_tmp+1:nar)
+                 endif
+ 
+              iw(1)=nar
+              do i=1,nar
+              iw(1+nar+i)=col(i)
+              enddo
+              if (nar > maxnar) stop 'increase sparsity fraction(spfra)'
+
+! CALLING IRLS TO SOLVE THE PROBLEM
+
+          leniw = 2*nar+1
+          lenrw = nar
+          dv = 0
+          atol = 1e-3
+          btol = 1e-3
+          conlim = 1200
+          itnlim = 1000
+          istop = 0
+          anorm = 0.0
+          acond = 0.0
+          arnorm = 0.0
+          xnorm = 0.0
+          localSize = n/4
+          
+        call LSMR(m, n, leniw, lenrw,iw,rw,cbst, damp,&
+           atol, btol, conlim, itnlim, localSize, nout,&
+           dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
+        if(istop==3) print*,'istop = 3, large condition number'
+
+
+        if (domain == 0.and.normt==0) then
+        do iiter = 1, initer-1
+        dtres=-cbst
+        call aprod(1,m,n,dv,dtres,leniw,lenrw,iw,rw)
+        do i=1,m
+        if(abs(dtres(i)).lt.tolr) then
+        wt(i)= 1.0/sqrt(abs(tolr))
+        else
+        wt(i)=1.0/sqrt(abs(dtres(i)))
+        endif
+        enddo
+        do i=1,nar
+        rw(i)=rw(i)*wt(iw(i+1))
+        enddo
+        do i=1,m
+        dtres(i)=cbst(i)*wt(i)
+        enddo
+
+          dv = 0
+          atol = 1e-3
+          btol = 1e-3
+          conlim = 1200
+          itnlim = 1000
+          istop = 0
+          anorm = 0.0
+          acond = 0.0
+          arnorm = 0.0
+          xnorm = 0.0
+          
+
+        call LSMR(m, n, leniw, lenrw,iw,rw,dtres, damp,&
+           atol, btol, conlim, itnlim, localSize, nout,&
+           dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
+        if(istop==3) print*,'istop = 3, large condition number'
+
+       do i=1,nar
+       rw(i)=rw(i)/wt(iw(i+1))
+       enddo
+
+       enddo ! finish inter interations for IRLS
+       endif
+    
+    
+    	if(domain==0.and.normt/=0) then
+        do iiter = 1, initer-1
+    	do i=1,n
+    	if (abs(dv(i)).lt.tolr) then
+    	rw(nar_tmp+i)=1.0/sqrt(tolr)*weight
+    	else
+    	rw(nar_tmp+i)=sqrt(1.0/abs(dv(i)))*weight
+    	endif
+    	enddo
+          dv = 0
+          atol = 1e-3
+          btol = 1e-3
+          conlim = 1200
+          itnlim = 1000
+          istop = 0
+          anorm = 0.0
+          acond = 0.0
+          arnorm = 0.0
+          xnorm = 0.0
+ 
+	   call LSMR(m, n, leniw, lenrw,iw,rw,cbst, damp,&
+           atol, btol, conlim, itnlim, localSize, nout,&
+           dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
+       if(istop==3) print*,'istop = 3, large condition number'
+          enddo
+     endif
+
+
+
+        if (domain/=0)then 
+          do iiter = 1,initer-1
+
+          dtres = 0
+          call aprod(1,m,n,dv,dtres,leniw,lenrw,iw,rw)
+          do i = nar_tmp+1,nar
+          if(abs(dtres(iw(1+i)))<tolr) then
+          rw(i)=1/sqrt(tolr)*rw(i+nars)
+          else
+          rw(i)=sqrt(1.0/abs(dtres(iw(1+i))))*rw(i+nars)
+          endif
+          enddo
+          dv = 0
+          atol = 1e-3
+          btol = 1e-3
+          conlim = 1200
+          itnlim = 1000
+          istop = 0
+          anorm = 0.0
+          acond = 0.0
+          arnorm = 0.0
+          xnorm = 0.0
+ 
+
+        call LSMR(m, n, leniw, lenrw,iw,rw,cbst, damp,&
+           atol, btol, conlim, itnlim, localSize, nout,&
+           dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
+        if(istop==3) print*,'istop = 3, large condition number'
+
+          enddo
+      endif
+
+       mean = sum(cbst(1:dall))/dall
+       std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
+       write(*,'(i2,a)'),iter,'th iteration...'
+       write(*,'(a,f7.3)'),'weight is:',weight
+       write(*,'(a,f8.1,a,f8.2,a,f8.3)'),'mean,std_devs and chi sqrue of &
+           residual: ',mean*1000,'ms ',1000*std_devs,'ms ',&
+           dnrm2(dall,cbst,1)**2/dall
+       write(66,'(i2,a)'),iter,'th iteration...'
+       write(66,'(a,f7.3)'),'weight is:',weight
+       write(66,'(a,f8.1,a,f8.2,a,f8.3)'),'mean,std_devs and chi sqrue of &
+           residual: ',mean*1000,'ms ',1000*std_devs,'ms ',&
+           dnrm2(dall,cbst,1)**2/dall
+
+        if (domain == 0) then
+    	call invwavetrans(nx-2,ny-2,nz-1,dv,maxlevel,maxleveld,&
+            HorizonType,VerticalType)
+        endif
+
+        write(*,'(a,2f7.4)'),'min and max velocity variation ',&
+            minval(dv),maxval(dv)
+        write(66,'(a,2f7.4)'),'min and max velocity variation ',&
+            minval(dv),maxval(dv)
+
+        do k=1,nz-1
+        do j=1,ny-2
+        do i=1,nx-2
+        if(dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i).ge.0.500) then
+        dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)=0.500
+        endif
+        if(dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i).le.-0.500) then
+        dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)=-0.500
+        endif
+        vsf(i+1,j+1,k)=vsf(i+1,j+1,k)+dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)
+       if(vsf(i+1,j+1,k).lt.Minvel) vsf(i+1,j+1,k)=Minvel
+       if(vsf(i+1,j+1,k).gt.Maxvel) vsf(i+1,j+1,k)=Maxvel
+        enddo
+        enddo
+        enddo
+        write(34,*)',OUTPUT S VELOCITY AT ITERATION',iter
+        do k=1,nz
+        do j=1,ny
+        write(34,'(100f7.3)') (vsf(i,j,k),i=1,nx)
+        enddo
+        enddo
+        write(34,*)',OUTPUT DWS AT ITERATION',iter
+        do k=1,nz-1
+        do j=2,ny-1
+        write(34,'(100f8.1)') (norm((k-1)*(ny-2)*(nx-2)+(j-2)*(nx-2)+i-1),i=2,nx-1)
+        enddo
+        enddo
+
+            enddo !end iteration
+
+! OUTPUT THE VELOCITY MODEL
+                
+        write(*,*),'Program finishes successfully'
+        write(66,*),'Program finishes successfully'
+
+        if(ifsyn == 1) then
+        open(65,file='Vs_model.real')
+	write(outsyn,'(a,a)') trim(inputfile),'Syn.dat'
+        open(63,file=outsyn)
+        do k=1,nz-1
+        do j=1,ny-2
+        do i=1,nx-2
+        write(65,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsftrue(i,j,k)
+        write(63,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i,j,k)
+        enddo
+        enddo
+        enddo
+        close(65)
+        close(63)
+        write(*,*),'Output True velocity model &
+                    to Vs_model.real'
+        write(*,*),'Output inverted shear velocity model &
+                    to ',outsyn
+        write(66,*),'Output True velocity model &
+                    to Vs_model.real'
+        write(66,*),'Output inverted shear velocity model &
+                    to ',outsyn
+        else
+	    write(outmodel,'(a,a)') trim(inputfile),'Measure.dat'
+        open(64,file=outmodel)
+        do k=1,nz-1
+        do j=1,ny-2
+        do i=1,nx-2
+        write(64,'(5f8.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k),vsf(i,j,k)
+        enddo
+        enddo
+        enddo
+        close(64)
+        write(*,*),'Output inverted shear velocity model &
+                    to ',outmodel
+        write(66,*),'Output inverted shear velocity model &
+                    to ',outmodel
+        endif
+
+        close(34)
+        close(40)
+        close(nout) !close lsmr.txt
+        close(66) !close surf_tomo.log
+        deallocate(obst)
+        deallocate(dsyn)
+        deallocate(dist)
+        deallocate(depz)
+        deallocate(scxf,sczf)
+        deallocate(rcxf,rczf)
+        deallocate(wavetype,igrt,nrc1)
+        deallocate(nsrc1,knum1,periods)
+        deallocate(rw)
+        deallocate(iw,col)
+        deallocate(cbst,wt,dtres,datweight)
+        deallocate(dv)
+        deallocate(norm)
+        deallocate(vsf)
+        deallocate(vsftrue)
+        if(kmaxRc.gt.0) then
+    	deallocate(tRc)
+    	endif
+    	if(kmaxRg.gt.0) then
+    	deallocate(tRg)
+    	endif
+    	if(kmaxLc.gt.0) then
+    	deallocate(tLc)
+    	endif
+    	if(kmaxLg.gt.0) then
+    	deallocate(tLg)
+    	endif
+
+            end program            

src/mainold.f90

@@ -0,0 +1,750 @@
+      ! CODE FOR SURFACE WAVE TOMOGRAPHY USING DISPERSION MEASUREMENTS
+      ! VERSION:
+      !      1.0 
+      ! AUTHOR:
+      !      HONGJIAN FANG. fanghj@mail.ustc.edu.cn
+      ! PURPOSE:
+      !      DIRECTLY INVERT SURFACE WAVE DISPERSION MEASUREMENTS FOR 3-D
+      ! STUCTURE WITHOUT THE INTERMEDIATE STEP OF CONSTUCTION THE PHASE
+      ! OR GROUP VELOCITY MAPS.
+      ! REFERENCE:
+      ! 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.
+      ! HISTORY:
+      !       2015/01/31 START TO REORGONIZE THE MESSY CODE   
+      !
+
+        program SurfTomo
+        use lsmrModule, only:lsmr
+	use  lsmrblasInterface, only : dnrm2
+        implicit none
+
+! VARIABLE DEFINE
+
+            character inputfile*80
+            character logfile*100
+            character outmodel*100
+            character outsyn*100
+            logical ex
+            character dummy*40
+            character datafile*80
+
+            integer nx,ny,nz
+            real goxd,gozd
+            real dvxd,dvzd
+            integer nsrc,nrc
+            real weight,weight0
+            real damp
+            real minthk
+            integer kmax,kmaxRc,kmaxRg,kmaxLc,kmaxLg
+            real*8,dimension(:),allocatable:: tRc,tRg,tLc,tLg
+            real,dimension(:),allocatable:: depz
+            integer itn
+            integer nout
+            integer localSize
+            real mean,std_devs,balances,balanceb
+            integer msurf
+            integer maxlevel,maxleveld
+            real,parameter:: tolr=1e-4
+            real,dimension(:),allocatable:: obst,dsyn,cbst,wt,dtres,dist,datweight
+         	real,dimension(:),allocatable:: pvall,depRp,pvRp
+        	real sta1_lat,sta1_lon,sta2_lat,sta2_lon
+        	real dist1
+                integer dall
+         	integer istep
+         	real,parameter :: pi=3.1415926535898
+         	integer checkstat
+         	integer ii,jj,kk
+                real, dimension (:,:), allocatable :: scxf,sczf
+                real, dimension (:,:,:), allocatable :: rcxf,rczf
+	        integer,dimension(:,:),allocatable::wavetype,igrt,nrc1
+         	integer,dimension(:),allocatable::nsrc1,knum1
+         	integer,dimension(:,:),allocatable::periods
+         	real,dimension(:),allocatable::rw
+         	integer,dimension(:),allocatable::iw,col
+                real,dimension(:),allocatable::dv,norm
+                real,dimension(:,:,:),allocatable::vsf
+                real,dimension(:,:,:),allocatable::vsftrue
+        	character strf
+        	integer veltp,wavetp
+        	real velvalue
+        	integer knum,knumo,err
+        	integer istep1,istep2
+        	integer period
+        	integer knumi,srcnum,count1
+        	integer HorizonType,VerticalType
+        	character line*200
+            integer iter,maxiter
+            integer iiter,initer
+            integer maxnar
+            real acond
+            real anorm
+            real arnorm
+            real rnorm
+            real xnorm
+            character str1
+            real atol,btol
+            real conlim
+            integer istop
+            integer itnlim
+            integer lenrw,leniw
+            integer nar,nar_tmp,nars
+            integer count3,nvz,nvx
+            integer m,maxvp,n
+            integer i,j,k
+            real Minvel,MaxVel
+            real spfra
+            integer domain,normt
+            real noiselevel
+            integer ifsyn
+            integer writepath
+            real averdws
+            real maxnorm
+	    real threshold,threshold0
+
+! OPEN FILES FIRST TO OUTPUT THE PROCESS
+            open(34,file='IterVel.out')
+            nout=36
+            open(nout,file='lsmr.txt')
+
+! OUTPUT PROGRAM INFOMATION            
+             write(*,*)
+             write(*,*),'                         S U R F  T O M O'
+             write(*,*),'PLEASE contact Hongjain Fang &
+                 (fanghj@mail.ustc.edu.cn) if you find any bug'
+             write(*,*)
+
+! READ INPUT FILE
+            if (iargc() < 1) then
+                write(*,*) 'input file [SurfTomo.in(default)]:'
+                read(*,'(a)') inputfile
+                if (len_trim(inputfile) <=1 ) then
+                    inputfile = 'SurfTomo.in'
+                else
+                    inputfile = inputfile(1:len_trim(inputfile))
+                endif
+            else
+                call getarg(1,inputfile)
+            endif
+            inquire(file = inputfile, exist = ex)
+            if (.not. ex) stop 'unable to open the inputfile'
+
+            open(10,file=inputfile,status='old')
+            read(10,'(a30)')dummy
+            read(10,'(a30)')dummy
+            read(10,'(a30)')dummy
+            read(10,*)datafile
+            read(10,*) nx,ny,nz
+            read(10,*) goxd,gozd
+            read(10,*) dvxd,dvzd
+            read(10,*) nsrc,nrc
+            read(10,*) weight0,damp
+            read(10,*) minthk
+            read(10,*) Minvel,Maxvel
+            read(10,*) domain,normt
+	        read(10,*) HorizonType,VerticalType
+            read(10,*) maxlevel,maxleveld
+            read(10,*) maxiter,initer
+            read(10,*) spfra
+            read(10,*) kmaxRc
+            write(*,*)  'model origin:latitude,longitue'
+            write(*,'(2f10.4)') goxd,gozd
+            write(*,*) 'grid spacing:latitude,longitue'
+            write(*,'(2f10.4)') dvxd,dvzd
+            write(*,*) 'model dimension:nx,ny,nz'
+            write(*,'(3i5)') nx,ny,nz
+	        write(logfile,'(a,a)')trim(inputfile),'.log' 
+            open(66,file=logfile)
+            write(66,*)
+            write(66,*),'                         S U R F  T O M O'
+            write(66,*),'PLEASE contact Hongjain Fang &
+                (fanghj@mail.ustc.edu.cn) if you find any bug'
+            write(66,*)
+            write(66,*) 'model origin:latitude,longitue'
+            write(66,'(2f10.4)') goxd,gozd
+            write(66,*) 'grid spacing:latitude,longitue'
+            write(66,'(2f10.4)') dvxd,dvzd
+            write(66,*) 'model dimension:nx,ny,nz'
+            write(66,'(3i5)') nx,ny,nz
+            if(kmaxRc.gt.0)then
+               allocate(tRc(kmaxRc),&
+                stat=checkstat)
+               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(66,*)'Rayleigh wave phase velocity used,periods:(s)'
+            write(66,'(50f6.2)')(tRc(i),i=1,kmaxRc)
+            endif
+            read(10,*)kmaxRg
+            if(kmaxRg.gt.0)then
+               allocate(tRg(kmaxRg), stat=checkstat)
+               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(66,*)'Rayleigh wave group velocity used,periods:(s)'
+            write(66,'(50f6.2)')(tRg(i),i=1,kmaxRg)
+            endif
+            read(10,*)kmaxLc
+            if(kmaxLc.gt.0)then
+               allocate(tLc(kmaxLc), stat=checkstat)
+               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(66,*)'Love wave phase velocity used,periods:(s)'
+            write(66,'(50f6.2)')(tLc(i),i=1,kmaxLc)
+            endif
+            read(10,*)kmaxLg
+            if(kmaxLg.gt.0)then
+               allocate(tLg(kmaxLg), stat=checkstat)
+               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(66,*)'Love wave group velocity used,periods:(s)'
+            write(66,'(50f6.2)')(tLg(i),i=1,kmaxLg)
+            endif
+            read(10,*)ifsyn
+            read(10,*)noiselevel
+	    read(10,*) threshold0
+            close(10)
+            kmax=kmaxRc+kmaxRg+kmaxLc+kmaxLg
+
+! READ MEASUREMENTS            
+            open(unit=87,file=datafile,status='old')
+            allocate(scxf(nsrc,kmax),sczf(nsrc,kmax),&
+            rcxf(nrc,nsrc,kmax),rczf(nrc,nsrc,kmax),stat=checkstat)
+            if(checkstat > 0)then
+            write(6,*)'error with allocate'
+            endif
+            allocate(periods(nsrc,kmax),wavetype(nsrc,kmax),&
+            nrc1(nsrc,kmax),nsrc1(kmax),knum1(kmax),&
+            igrt(nsrc,kmax),stat=checkstat)
+            if(checkstat > 0)then
+            write(6,*)'error with allocate'
+            endif
+            allocate(obst(nrc*nsrc*kmax),dist(nrc*nsrc*kmax),&
+            stat=checkstat)
+            allocate(pvall(nrc*nsrc*kmax),depRp(nrc*nsrc*kmax),&
+            pvRp(nrc*nsrc*kmax),stat=checkstat)
+            IF(checkstat > 0)THEN
+            write(6,*)'error with allocate'
+            ENDIF
+            istep=0
+            istep2=0
+            dall=0
+            knumo=12345
+	        knum=0
+	        istep1=0
+            do 
+            read(87,'(a)',iostat=err) line
+            if(err.eq.0) then
+            if(line(1:1).eq.'#') then
+            read(line,*) str1,sta1_lat,sta1_lon,period,wavetp,veltp
+            if(wavetp.eq.2.and.veltp.eq.0) knum=period
+            if(wavetp.eq.2.and.veltp.eq.1) knum=kmaxRc+period
+            if(wavetp.eq.1.and.veltp.eq.0) knum=kmaxRg+kmaxRc+period
+            if(wavetp.eq.1.and.veltp.eq.1) knum=kmaxLc+kmaxRg+&
+               kmaxRc+period
+            if(knum.ne.knumo) then
+            istep=0
+            istep2=istep2+1
+            endif
+            istep=istep+1
+            istep1=0
+            sta1_lat=(90.0-sta1_lat)*pi/180.0
+            sta1_lon=sta1_lon*pi/180.0
+            scxf(istep,knum)=sta1_lat
+            sczf(istep,knum)=sta1_lon
+            periods(istep,knum)=period
+            wavetype(istep,knum)=wavetp
+            igrt(istep,knum)=veltp
+            nsrc1(knum)=istep
+            knum1(istep2)=knum
+            knumo=knum
+            else
+            read(line,*) sta2_lat,sta2_lon,velvalue
+            istep1=istep1+1
+            dall=dall+1
+            sta2_lat=(90.0-sta2_lat)*pi/180.0
+            sta2_lon=sta2_lon*pi/180.0
+            rcxf(istep1,istep,knum)=sta2_lat
+            rczf(istep1,istep,knum)=sta2_lon
+            call delsph(sta1_lat,sta1_lon,sta2_lat,sta2_lon,dist1)
+	    dist(dall)=dist1
+            obst(dall)=dist1/velvalue
+            pvall(dall)=velvalue
+            nrc1(istep,knum)=istep1
+            endif
+            else
+            exit
+            endif
+            enddo
+            close(87)
+            allocate(depz(nz), stat=checkstat)
+            maxnar = dall*nx*ny*nz*spfra!sparsity fraction
+            maxvp = (nx-2)*(ny-2)*(nz-1)
+            allocate(dv(maxvp), stat=checkstat)
+            allocate(norm(maxvp), stat=checkstat)
+            allocate(vsf(nx,ny,nz), stat=checkstat)
+            allocate(vsftrue(nx,ny,nz), stat=checkstat)
+        	
+        	allocate(rw(maxnar), stat=checkstat)
+        	if(checkstat > 0)then
+        	   write(6,*)'error with allocate: real rw'
+        	endif
+        	allocate(iw(2*maxnar+1), stat=checkstat)
+        	if(checkstat > 0)then
+        	   write(6,*)'error with allocate: integer iw'
+        	endif
+        	allocate(col(maxnar), stat=checkstat)
+        	if(checkstat > 0)then
+        	   write(6,*)'error with allocate:  integer iw'
+        	endif
+           allocate(cbst(dall+maxvp),dsyn(dall),datweight(dall),wt(dall+maxvp),dtres(dall+maxvp),&
+                stat=checkstat)
+
+! MEASUREMENTS STATISTICS AND READ INITIAL MODEL           
+            write(*,'(a,i7)') 'Number of all measurements',dall
+
+            open(10,file='MOD',status='old')
+            read(10,*) (depz(i),i=1,nz)
+            do k = 1,nz
+            do j = 1,ny
+            read(10,*)(vsf(i,j,k),i=1,nx)
+            enddo
+            enddo
+            close(10)
+            write(*,*) 'grid points in depth direction:(km)'
+            write(*,'(50f6.2)') depz
+
+
+
+! CHECKERBOARD TEST
+            if (ifsyn == 1) then
+            write(*,*) 'Checkerboard Resolution Test Begin'
+            vsftrue = vsf
+
+            open(11,file='MOD.true',status='old')
+            do k = 1,nz
+            do j = 1,ny
+            read(11,*) (vsftrue(i,j,k),i=1,nx)
+            enddo
+            enddo
+            close(11)
+
+            call synthetic(nx,ny,nz,maxvp,vsftrue,obst,&
+                goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
+                tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
+                scxf,sczf,rcxf,rczf,nrc1,nsrc1,knum1,kmax,&
+                nsrc,nrc,noiselevel)
+            endif
+
+
+
+! ITERATE UNTILL CONVERGE
+            writepath = 0
+            do iter = 1,maxiter
+            iw = 0
+            rw = 0.0
+            col = 0
+
+! COMPUTE SENSITIVITY MATRIX
+            if (iter == maxiter) then
+                writepath = 1
+                open(40,file='raypath.out')
+            endif
+            write(*,*) 'computing sensitivity matrix...'
+            call CalSurfG(nx,ny,nz,maxvp,vsf,iw,rw,col,dsyn,&
+                goxd,gozd,dvxd,dvzd,kmaxRc,kmaxRg,kmaxLc,kmaxLg,&
+                tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk,&
+                scxf,sczf,rcxf,rczf,nrc1,nsrc1,knum1,kmax,&
+                nsrc,nrc,nar,domain,&
+                maxlevel,maxleveld,HorizonType,VerticalType,writepath)
+
+	        do i = 1,dall
+	        cbst(i) = obst(i) - dsyn(i)
+	        enddo
+
+		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
+		cbst(i) = cbst(i)*datweight(i)
+		enddo
+
+		do i = 1,nar
+		rw(i) = rw(i)*datweight(iw(1+i))
+		enddo
+ 
+            	norm=0
+            	do i=1,nar
+            	norm(col(i))=norm(col(i))+abs(rw(i))
+            	enddo
+            	averdws=0
+            	maxnorm=0
+            	do i=1,maxvp
+            	averdws =  averdws+norm(i)
+            	if(norm(i)>maxnorm) maxnorm=norm(i)
+            	enddo
+            	averdws=averdws/maxvp
+            	write(66,*)'Maximum and Average DWS values:',maxnorm,averdws
+            	write(66,*)'Threshold is:',threshold
+            
+! WRITE OUT RESIDUAL FOR THE FIRST AND LAST ITERATION
+               if(iter.eq.1) then
+               open(88,file='residualFirst.dat')
+               do i=1,dall
+               write(88,*) dist(i),dsyn(i),obst(i), &
+               dsyn(i)*datweight(i),obst(i)*datweight(i),datweight(i)
+               enddo
+               close(88)
+               endif
+               if(iter.eq.maxiter) then
+               open(88,file='residualLast.dat')
+               do i=1,dall
+               write(88,*) dist(i),dsyn(i),obst(i), &
+               dsyn(i)*datweight(i),obst(i)*datweight(i),datweight(i)
+               enddo
+               close(88)
+               endif
+ 
+         
+! ADDING REGULARIZATION TERM
+       	       weight=dnrm2(dall,cbst,1)**2/dall*weight0
+               nar_tmp=nar
+     	       nars=0
+               if (domain == 0 .and. normt==0) then
+               do i=1,maxvp
+               rw(nar+i)=weight
+               iw(1+nar+i)=dall+i
+               col(nar+i)=i
+               cbst(dall+i)=0
+               enddo
+               nar = nar + maxvp
+                m = dall + maxvp
+                n = maxvp
+             	elseif(domain == 0 .and. normt/=0) then
+             	do i=1,maxvp
+                rw(nar+i)=weight
+                iw(1+nar+i)=dall+i
+                col(nar+i)=i
+                cbst(dall+i)=0
+                enddo
+                nar = nar + maxvp
+                 m = dall + maxvp
+                 n = maxvp
+                 else
+                 count3=0
+                 nvz=ny-2
+                 nvx=nx-2
+                 do k=1,nz-1
+                 do j=1,nvz
+                 do i=1,nvx
+                 if(i==1.or.i==nvx.or.j==1.or.j==nvz.or.k==1.or.k==nz-1)then
+                 count3=count3+1
+                 col(nar+1)=(k-1)*nvz*nvx+(j-1)*nvx+i
+                 rw(nar+1)=2.0*weight
+                 iw(1+nar+1)=dall+count3
+                 cbst(dall+count3)=0
+                 nar=nar+1
+                 else
+                 count3=count3+1
+                 col(nar+1)=(k-1)*nvz*nvx+(j-1)*nvx+i
+                 rw(nar+1)=6.0*weight
+                 iw(1+nar+1)=dall+count3
+                 rw(nar+2)=-1.0*weight
+                 iw(1+nar+2)=dall+count3
+                 col(nar+2)=(k-1)*nvz*nvx+(j-1)*nvx+i-1
+                 rw(nar+3)=-1.0*weight
+                 iw(1+nar+3)=dall+count3
+                 col(nar+3)=(k-1)*nvz*nvx+(j-1)*nvx+i+1
+                 rw(nar+4)=-1.0*weight
+                 iw(1+nar+4)=dall+count3
+                 col(nar+4)=(k-1)*nvz*nvx+(j-2)*nvx+i
+                 rw(nar+5)=-1.0*weight
+                 iw(1+nar+5)=dall+count3
+                 col(nar+5)=(k-1)*nvz*nvx+j*nvx+i
+                 rw(nar+6)=-1.0*weight
+                 iw(1+nar+6)=dall+count3
+                 col(nar+6)=(k-2)*nvz*nvx+(j-1)*nvx+i
+                 rw(nar+7)=-1.0*weight
+                 iw(1+nar+7)=dall+count3
+                 col(nar+7)=k*nvz*nvx+(j-1)*nvx+i
+                 cbst(dall+count3)=0
+                 nar=nar+7
+                 endif
+                 enddo
+                 enddo
+                 enddo
+                   m = dall + count3
+                   n = maxvp
+                   nars = nar - nar_tmp
+                   rw(nar+1:nar+nars) = rw(nar_tmp+1:nar)
+                 endif
+ 
+              iw(1)=nar
+              do i=1,nar
+              iw(1+nar+i)=col(i)
+              enddo
+              if (nar > maxnar) stop 'increase sparsity fraction(spfra)'
+
+! CALLING IRLS TO SOLVE THE PROBLEM
+
+          leniw = 2*nar+1
+          lenrw = nar
+          dv = 0
+          atol = 1e-3
+          btol = 1e-3
+          conlim = 1200
+          itnlim = 1000
+          istop = 0
+          anorm = 0.0
+          acond = 0.0
+          arnorm = 0.0
+          xnorm = 0.0
+          localSize = n/4
+          
+        call LSMR(m, n, leniw, lenrw,iw,rw,cbst, damp,&
+           atol, btol, conlim, itnlim, localSize, nout,&
+           dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
+        if(istop==3) print*,'istop = 3, large condition number'
+
+
+        if (domain == 0.and.normt==0) then
+        do iiter = 1, initer
+        dtres=-cbst
+        call aprod(1,m,n,dv,dtres,leniw,lenrw,iw,rw)
+        do i=1,m
+        if(abs(dtres(i)).lt.tolr) then
+        wt(i)= 1.0/sqrt(abs(tolr))
+        else
+        wt(i)=1.0/sqrt(abs(dtres(i)))
+        endif
+        enddo
+        do i=1,nar
+        rw(i)=rw(i)*wt(iw(i+1))
+        enddo
+        do i=1,m
+        dtres(i)=cbst(i)*wt(i)
+        enddo
+
+          dv = 0
+          atol = 1e-3
+          btol = 1e-3
+          conlim = 1200
+          itnlim = 1000
+          istop = 0
+          anorm = 0.0
+          acond = 0.0
+          arnorm = 0.0
+          xnorm = 0.0
+          
+
+        call LSMR(m, n, leniw, lenrw,iw,rw,dtres, damp,&
+           atol, btol, conlim, itnlim, localSize, nout,&
+           dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
+        if(istop==3) print*,'istop = 3, large condition number'
+
+       do i=1,nar
+       rw(i)=rw(i)/wt(iw(i+1))
+       enddo
+
+       enddo ! finish inter interations for IRLS
+       endif
+    
+    
+    	if(domain==0.and.normt/=0) then
+        do iiter = 1, initer
+    	do i=1,n
+    	if (abs(dv(i)).lt.tolr) then
+    	rw(nar_tmp+i)=1.0/sqrt(tolr)*weight
+    	else
+    	rw(nar_tmp+i)=sqrt(1.0/abs(dv(i)))*weight
+    	endif
+    	enddo
+          dv = 0
+          atol = 1e-3
+          btol = 1e-3
+          conlim = 1200
+          itnlim = 1000
+          istop = 0
+          anorm = 0.0
+          acond = 0.0
+          arnorm = 0.0
+          xnorm = 0.0
+ 
+	   call LSMR(m, n, leniw, lenrw,iw,rw,cbst, damp,&
+           atol, btol, conlim, itnlim, localSize, nout,&
+           dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
+       if(istop==3) print*,'istop = 3, large condition number'
+          enddo
+     endif
+
+
+
+        if (domain/=0)then 
+          do iiter = 1,initer
+
+          dtres = 0
+          call aprod(1,m,n,dv,dtres,leniw,lenrw,iw,rw)
+          do i = nar_tmp+1,nar
+          if(abs(dtres(iw(1+i)))<tolr) then
+          rw(i)=1/sqrt(tolr)*rw(i+nars)
+          else
+          rw(i)=sqrt(1.0/abs(dtres(iw(1+i))))*rw(i+nars)
+          endif
+          enddo
+          dv = 0
+          atol = 1e-3
+          btol = 1e-3
+          conlim = 1200
+          itnlim = 1000
+          istop = 0
+          anorm = 0.0
+          acond = 0.0
+          arnorm = 0.0
+          xnorm = 0.0
+ 
+
+        call LSMR(m, n, leniw, lenrw,iw,rw,cbst, damp,&
+           atol, btol, conlim, itnlim, localSize, nout,&
+           dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
+        if(istop==3) print*,'istop = 3, large condition number'
+
+          enddo
+      endif
+
+       mean = sum(cbst(1:dall))/dall
+       std_devs = sqrt(sum(cbst(1:dall)**2)/dall - mean**2)
+       write(*,'(i2,a)'),iter,'th iteration...'
+       write(*,'(a,f7.3)'),'weight is:',weight
+       write(*,'(a,f8.1,a,f8.2,a,f8.3)'),'mean,std_devs and chi sqrue of &
+           residual: ',mean*1000,'ms ',1000*std_devs,'ms ',&
+           dnrm2(dall,cbst,1)**2/dall
+       write(66,'(i2,a)'),iter,'th iteration...'
+       write(66,'(a,f7.3)'),'weight is:',weight
+       write(66,'(a,f8.1,a,f8.2,a,f8.3)'),'mean,std_devs and chi sqrue of &
+           residual: ',mean*1000,'ms ',1000*std_devs,'ms ',&
+           dnrm2(dall,cbst,1)**2/dall
+
+        if (domain == 0) then
+    	call invwavetrans(nx-2,ny-2,nz-1,dv,maxlevel,maxleveld,&
+            HorizonType,VerticalType)
+        endif
+
+        write(*,'(a,2f7.4)'),'min and max velocity variation ',&
+            minval(dv),maxval(dv)
+        write(66,'(a,2f7.4)'),'min and max velocity variation ',&
+            minval(dv),maxval(dv)
+
+        do k=1,nz-1
+        do j=1,ny-2
+        do i=1,nx-2
+        if(dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i).ge.0.500) then
+        dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)=0.500
+        endif
+        if(dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i).le.-0.500) then
+        dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)=-0.500
+        endif
+        vsf(i+1,j+1,k)=vsf(i+1,j+1,k)+dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)
+       if(vsf(i+1,j+1,k).lt.Minvel) vsf(i+1,j+1,k)=Minvel
+       if(vsf(i+1,j+1,k).gt.Maxvel) vsf(i+1,j+1,k)=Maxvel
+        enddo
+        enddo
+        enddo
+        write(34,*)',OUTPUT S VELOCITY AT ITERATION',iter
+        do k=1,nz
+        do j=1,ny
+        write(34,'(100f7.3)') (vsf(i,j,k),i=1,nx)
+        enddo
+        enddo
+        write(34,*)',OUTPUT DWS AT ITERATION',iter
+        do k=1,nz-1
+        do j=2,ny-1
+        write(34,'(100f8.1)') (norm((k-1)*(ny-2)*(nx-2)+(j-2)*(nx-2)+i-1),i=2,nx-1)
+        enddo
+        enddo
+
+            enddo !end iteration
+
+! OUTPUT THE VELOCITY MODEL
+                
+        write(*,*),'Program finishes successfully'
+        write(66,*),'Program finishes successfully'
+
+        if(ifsyn == 1) then
+        open(65,file='Vs_model.real')
+	    write(outsyn,'(a,a)') trim(inputfile),'Syn.dat'
+        open(63,file=outsyn)
+        do k=1,nz
+        do j=1,ny
+        write(65,'(100f7.3)') (vsftrue(i,j,k),i=1,nx)
+        write(63,'(100f7.3)') (vsf(i,j,k),i=1,nx)
+        enddo
+        enddo
+        close(65)
+        close(63)
+        write(*,*),'Output True velocity model &
+                    to Vs_model.real'
+        write(*,*),'Output inverted shear velocity model &
+                    to ',outsyn
+        write(66,*),'Output True velocity model &
+                    to Vs_model.real'
+        write(66,*),'Output inverted shear velocity model &
+                    to ',outsyn
+        else
+	    write(outmodel,'(a,a)') trim(inputfile),'Measure.dat'
+        open(64,file=outmodel)
+        do k=1,nz
+        do j=1,ny
+        write(64,'(100f7.3)') (vsf(i,j,k),i=1,nx)
+        enddo
+        enddo
+        close(64)
+        write(*,*),'Output inverted shear velocity model &
+                    to ',outmodel
+        write(66,*),'Output inverted shear velocity model &
+                    to ',outmodel
+        endif
+
+        close(34)
+        close(40)
+        close(nout) !close lsmr.txt
+        close(66) !close surf_tomo.log
+        deallocate(obst)
+        deallocate(dsyn)
+        deallocate(dist)
+        deallocate(depz)
+        deallocate(scxf,sczf)
+        deallocate(rcxf,rczf)
+        deallocate(wavetype,igrt,nrc1)
+        deallocate(nsrc1,knum1,periods)
+        deallocate(rw)
+        deallocate(iw,col)
+        deallocate(cbst,wt,dtres,datweight)
+        deallocate(dv)
+        deallocate(norm)
+        deallocate(vsf)
+        deallocate(vsftrue)
+        if(kmaxRc.gt.0) then
+    	deallocate(tRc)
+    	endif
+    	if(kmaxRg.gt.0) then
+    	deallocate(tRg)
+    	endif
+    	if(kmaxLc.gt.0) then
+    	deallocate(tLc)
+    	endif
+    	if(kmaxLg.gt.0) then
+    	deallocate(tLg)
+    	endif
+
+            end program            

src/merge1.f90

@@ -0,0 +1,21 @@
+	subroutine merge1(vec,n)
+	implicit none
+	integer n
+	real vec(n)
+	
+	integer half,start,endt
+	integer i
+	real tmp
+	half=n/2
+   	start=half
+	endt=half+1
+	do while(start.gt.1)
+	do i=start,endt-1,2
+	tmp=vec(i)
+	vec(i)=vec(i+1)
+	vec(i+1)=tmp
+	enddo
+	start=start-1
+	endt=endt+1
+	enddo
+	end subroutine

src/split.f90

@@ -0,0 +1,24 @@
+	subroutine split(vec,n)
+	implicit none
+	integer n
+	real vec(n)
+	
+!c local 
+	integer start,endt,i,j
+	real tmp
+	start=2
+	!if(mod(n,2).ne.0) then
+	! endt=n-1
+	!else
+	endt=n
+	!endif
+	do while (start.lt.endt)
+	do i=start,endt-1,2
+	tmp=vec(i)
+	vec(i)=vec(i+1)
+	vec(i+1)=tmp
+	end do
+	start=start+1
+	endt=endt-1
+	enddo
+	end subroutine

src/waveletD8.f90

@@ -0,0 +1,113 @@
+        subroutine transformD8(a,n)
+        implicit none
+        integer n
+        real a(n)
+
+        integer i,j
+        integer half
+        real tmp(n)
+        real*8 h(8),g(8)
+        data h/-0.010597401784997278,&
+        0.032883011666982945,&
+        0.030841381835986965,&
+        -0.18703481171888114,&
+        -0.02798376941698385,&
+        0.6308807679295904,&
+        0.7148465705525415,&
+        0.23037781330885523/
+
+
+        data g/-0.23037781330885523,&
+        0.7148465705525415,&
+        -0.6308807679295904,&
+        -0.02798376941698385,&
+        0.18703481171888114,&
+        0.030841381835986965,&
+        -0.032883011666982945,&
+        -0.010597401784997278/
+
+        half=n/2
+        i=1
+        tmp=0
+        do j=1,n-7,2
+        tmp(i)=a(j)*h(1)+a(j+1)*h(2)+a(j+2)*h(3)+a(j+3)*h(4)+a(j+4)*h(5) &
+        +a(j+5)*h(6)+a(j+6)*h(7)+a(j+7)*h(8)
+        tmp(i+half)=a(j)*g(1)+a(j+1)*g(2)+a(j+2)*g(3)+a(j+3)*g(4) &
+        +a(j+4)*g(5)+a(j+5)*g(6)+a(j+6)*g(7)+a(j+7)*g(8)
+        i=i+1
+        enddo
+
+        tmp(i)=a(n-5)*h(1)+a(n-4)*h(2)+a(n-3)*h(3)+a(n-2)*h(4)+a(n-1)*h(5) &
+        +a(n)*h(6)+a(1)*h(7)+a(2)*h(8)
+        tmp(i+half)=a(n-5)*g(1)+a(n-4)*g(2)+a(n-3)*g(3)+a(n-2)*g(4) &
+        +a(n-1)*g(5) +a(n)*g(6)+a(1)*g(7)+a(2)*g(8)
+        tmp(i+1)=a(n-3)*h(1)+a(n-2)*h(2)+a(n-1)*h(3)+a(n)*h(4)+a(1)*h(5) &
+        +a(2)*h(6)+a(3)*h(7)+a(4)*h(8)
+        tmp(i+1+half)=a(n-3)*g(1)+a(n-2)*g(2)+a(n-1)*g(3)+a(n)*g(4) &
+        +a(1)*g(5)+a(2)*g(6)+a(3)*g(7)+a(4)*g(8)
+        tmp(i+2)=a(n-1)*h(1)+a(n)*h(2)+a(1)*h(3)+a(2)*h(4)+a(3)*h(5) &
+        +a(4)*h(6)+a(5)*h(7)+a(6)*h(8)
+        tmp(i+2+half)=a(n-1)*g(1)+a(n)*g(2)+a(1)*g(3)+a(2)*g(4)+a(3)*g(5) &
+        +a(4)*g(6)+a(5)*g(7)+a(6)*g(8)
+
+        do i=1,n
+        a(i)=tmp(i)
+        enddo
+         
+        end subroutine
+
+        subroutine invTransformD8(a,n)
+        implicit none
+        integer n
+        real a(n)
+        
+        real tmp(n)
+        real*8 Ih(8),Ig(8)
+        integer half
+        integer i,j
+        data Ih/0.23037781330885523,&
+        0.7148465705525415,&
+         0.6308807679295904,&
+        -0.02798376941698385,&
+        -0.18703481171888114,&
+        0.030841381835986965,&
+        0.032883011666982945,&
+        -0.010597401784997278/
+        data Ig/-0.010597401784997278,&
+        -0.032883011666982945,&
+        0.030841381835986965,&
+        0.18703481171888114,&
+        -0.02798376941698385,&
+        -0.6308807679295904,&
+        0.7148465705525415,&
+        -0.23037781330885523/
+
+        half=n/2
+        tmp(2)=a(half-2)*Ih(1)+a(n-2)*Ig(1)+a(half-1)*Ih(3)+a(n-1)*Ig(3) &
+        +a(half)*Ih(5)+a(n)*Ig(5)+a(1)*Ih(7)+a(half+1)*Ig(7)
+        tmp(1)=a(half-2)*Ih(2)+a(n-2)*Ig(2)+a(half-1)*Ih(4)+a(n-1)*Ig(4) &
+        +a(half)*Ih(6)+a(n)*Ig(6)+a(1)*Ih(8)+a(half+1)*Ig(8)
+        tmp(4)=a(half-1)*Ih(1)+a(n-1)*Ig(1)+a(half)*Ih(3)+a(n)*Ig(3) &
+        +a(1)*Ih(5)+a(half+1)*Ig(5)+a(2)*Ih(7)+a(half+2)*Ig(7)
+        tmp(3)=a(half-1)*Ih(2)+a(n-1)*Ig(2)+a(half)*Ih(4)+a(n)*Ig(4) &
+        +a(1)*Ih(6)+a(half+1)*Ig(6)+a(2)*Ih(8)+a(half+2)*Ig(8)
+        tmp(6)=a(half)*Ih(1)+a(n)*Ig(1)+a(1)*Ih(3)+a(half+1)*Ig(3) &
+        +a(2)*Ih(5)+a(half+2)*Ig(5)+a(3)*Ih(7)+a(half+3)*Ig(7)
+        tmp(5)=a(half)*Ih(2)+a(n)*Ig(2)+a(1)*Ih(4)+a(half+1)*Ig(4) &
+        +a(2)*Ih(6)+a(half+2)*Ig(6)+a(3)*Ih(8)+a(half+3)*Ig(8)
+
+        j=6
+        do i=1,half-3
+        j=j+1
+        tmp(j)=a(i)*Ih(2)+a(i+half)*Ig(2)+a(i+1)*Ih(4)+a(i+1+half)*Ig(4) &
+        +a(i+2)*Ih(6)+a(i+2+half)*Ig(6)+a(i+3)*Ih(8)+a(i+3+half)*Ig(8)
+        j=j+1
+        tmp(j)=a(i)*Ih(1)+a(i+half)*Ig(1)+a(i+1)*Ih(3)+a(i+1+half)*Ig(3) &
+        +a(i+2)*Ih(5)+a(i+2+half)*Ig(5)+a(i+3)*Ih(7)+a(i+3+half)*Ig(7)
+        enddo
+
+        do i=1,n
+        a(i)=tmp(i)
+        enddo
+
+        end subroutine

src/wavelettrans3domp.f90

@@ -0,0 +1,90 @@
+	subroutine wavelettrans(nx,ny,nz,row,maxlevel,maxleveld,HorizonType,VerticalType)
+        use omp_lib
+	implicit none
+	integer nx,ny,nz,maxlevel,maxleveld
+	real row(*)
+
+	integer j,k
+	real fxs(nx),fzs(nz),fys(ny)
+	integer HorizonType,VerticalType
+!	if(PorS == 1 .or. PorS == 3) then
+!!$omp parallel & 
+!!$omp default(private) &
+!!$omp shared(nz,ny,nx,maxlevel,row,HorizonType) 
+!!$omp do
+!	do k=1,nz
+!	do j=1,ny
+!	fxp=row(1+(j-1)*nx+(k-1)*nx*ny:nx+(j-1)*nx+(k-1)*nx*ny)
+!	call forwardtrans(fxp,nx,maxlevel,HorizonType)
+!	row(1+(j-1)*nx+(k-1)*nx*ny:nx+(j-1)*nx+(k-1)*nx*ny)=fxp
+!	enddo
+!	enddo
+!!$omp end do
+!!$omp end parallel
+!!$omp parallel & 
+!!$omp default(private) &
+!!$omp shared(nz,ny,nx,maxlevel,row,HorizonType) 
+!!$omp do
+!	do k=1,nz
+!	do j=1,nx
+!	fyp=row(j+(k-1)*nx*ny:j+(ny-1)*nx+(k-1)*nx*ny:nx)
+!	call forwardtrans(fyp,ny,maxlevel,HorizonType)
+!	row(j+(k-1)*nx*ny:j+(ny-1)*nx+(k-1)*nx*ny:nx)=fyp
+!	enddo
+!	enddo
+!!$omp end do
+!!$omp end parallel
+!!$omp parallel & 
+!!$omp default(private) &
+!!$omp shared(nz,ny,nx,maxleveld,row,VerticalType) 
+!!$omp do
+!	do k=1,ny
+!	do j=1,nx
+!        fzp=row(j+(k-1)*nx:j+(k-1)*nx+(nz-1)*nx*ny:nx*ny)
+!        call forwardtrans(fzp,nz,maxleveld,VerticalType)
+!        row(j+(k-1)*nx:j+(k-1)*nx+(nz-1)*nx*ny:nx*ny)=fzp
+!	enddo
+!	enddo
+!!$omp end do
+!!$omp end parallel
+!	endif
+!$omp parallel & 
+!$omp default(private) &
+!$omp shared(nz,ny,nx,maxlevel,row,HorizonType) 
+!$omp do
+	do k=1,nz
+	do j=1,ny
+	fxs=row(1+(j-1)*nx+(k-1)*nx*ny:nx+(j-1)*nx+(k-1)*nx*ny)
+	call forwardtrans(fxs,nx,maxlevel,HorizonType)
+	row(1+(j-1)*nx+(k-1)*nx*ny:nx+(j-1)*nx+(k-1)*nx*ny)=fxs
+	enddo
+	enddo
+!$omp end do
+!$omp end parallel
+!$omp parallel & 
+!$omp default(private) &
+!$omp shared(nz,ny,nx,maxlevel,row,HorizonType) 
+!$omp do
+	do k=1,nz
+	do j=1,nx
+	fys=row(j+(k-1)*nx*ny:j+(ny-1)*nx+(k-1)*nx*ny:nx)  
+	call forwardtrans(fys,ny,maxlevel,HorizonType)
+	row(j+(k-1)*nx*ny:j+(ny-1)*nx+(k-1)*nx*ny:nx)=fys
+	enddo
+	enddo
+!$omp end do
+!$omp end parallel
+!$omp parallel & 
+!$omp default(private) &
+!$omp shared(nz,ny,nx,maxleveld,row,VerticalType) 
+!$omp do
+	 do k=1,ny
+	  do j=1,nx
+		fzs=row(j+(k-1)*nx:j+(k-1)*nx+(nz-1)*nx*ny:nx*ny)
+		call forwardtrans(fzs,nz,maxleveld,VerticalType)
+		row(j+(k-1)*nx:j+(k-1)*nx+(nz-1)*nx*ny:nx*ny)=fzs
+	  enddo
+	 enddo
+!$omp end do
+!$omp end parallel
+	end subroutine