CODING THE TOMOGRAPHY

I am using subroutine xtomo() .

subroutine xtomo( adj,add,z0,dz,x0,dx,h0,dh,y0,dy,modl,nt,nz,nx,  data,   nh,ny)
integer           adj,add,                             nt,nz,nx,          nh,ny
real		                                  modl(nt,nz,nx), data(nt,nh,ny)
integer	                                               it,iz,ix,          ih,iy
real                      z0,dz,x0,dx,h0,dh,y0,dy
real		          z,    x,    h,    y,                zmax
call adjnull( adj,add, modl, nt*nz*nx, data, nt*nh*ny);       zmax=z0+dz*(nz-1)
do ix= 1, nx {  x= x0 + dx*(ix-1)
do ih= 1, nh {	h= h0 + dh*(ih-1)
do iz= 1, nz {	z= z0 + dz*(iz-1)

y = x +  h *  z/zmax    # maybe opposite sign on h or (1-z/zmax)

iy = 1.5 + (y-y0) / dy
	if( 1 <= iy & iy <= ny)
		do it= 1, nt
			if( adj == 0)
				data(it,ih,iy) = data(it,ih,iy) + modl(it,iz,ix)
			else
				modl(it,iz,ix) = modl(it,iz,ix) + data(it,ih,iy)
	}}}
return; end

The inversion problem in tomography was solved by using the conjugate-gradient subroutine cgstep() described in PVI invoked by subroutine tominv() . It allows for damping eps but I did not test damping.

subroutine tominv( niter,modl,nt,nz,nx,  data,   nh,ny, z0,dz,x0,dx,h0,dh,y0,dy)
integer      iter, niter,     nt,nz,nx,          nh,ny
	real	         modl(nt,nz,nx), data(nt,nh,ny), eps
	real	                 z0,x0,          h0,y0
	real	                 dz,dx,          dh,dy
integer                nmodl,           ndata,            ir
temporary real                             rr(nt*nh*ny + nt*nz*nx)
temporary real          dmodl(nt,nz,nx),   dr(nt*nh*ny + nt*nz*nx)
temporary real          smodl(nt,nz,nx),   sr(nt*nh*ny + nt*nz*nx)
                       nmodl= nt*nz*nx; ndata=nt*nh*ny;  ir= ndata + 1
call null(  modl, nmodl )
call copy(            ndata, data, rr( 1));	  eps= 0.  # Turn off.
call ident( 0,0, eps, nmodl, modl, rr(ir))
do iter= 0, niter {
	call snap( 'rr.H',   nt*nh, ny, rr)			# snapshot
	call xtomo( 1,0, z0,dz,x0,dx,h0,dh,y0,dy, dmodl, nt,nz,nx,rr( 1), nh,ny)
	call ident( 1,1, eps, nmodl,              dmodl,          rr(ir)       )
	call xtomo( 0,0, z0,dz,x0,dx,h0,dh,y0,dy, dmodl, nt,nz,nx,dr( 1), nh,ny)
	call ident( 0,0, eps, nmodl,              dmodl,          dr(ir)       )
	call cgstep( iter,       nmodl, modl, dmodl,smodl,
			   ndata+nmodl,   rr, dr,   sr    )
	call snap( 'modl.H', nt*nz, nx, modl)			# snapshot
	}
return; end

I omitted the listing of the main program because it merely serves to connect the listed subroutines to our file system.