diff --git a/bin/DSurfTomo b/bin/DSurfTomo deleted file mode 100755 index 00a895b..0000000 Binary files a/bin/DSurfTomo and /dev/null differ diff --git a/bin/DSurfTomo b/bin/DSurfTomo new file mode 120000 index 0000000..0b35a8a --- /dev/null +++ b/bin/DSurfTomo @@ -0,0 +1 @@ +../src/DSurfTomo \ No newline at end of file diff --git a/src/CalSurfG.f90 b/src/CalSurfG.f90 index 06b1005..470ccc4 100644 --- a/src/CalSurfG.f90 +++ b/src/CalSurfG.f90 @@ -92,6 +92,17 @@ subroutine depthkernel(nx,ny,nz,vel,pvRc,sen_vsRc,sen_vpRc,sen_rhoRc, & dlncg_dlnvs(nn,i) = (cg2(nn)-cg1(nn))/(dlnVs*vsz(i)) enddo + !note here, no integral from z1 to z2 is needed. The grid based nodes + ! have already taken into consideration of it. Layer based need integration. + ! Test passed + !if (i == 1) then + ! dlncg_dlnvs(1:kmaxRc,i) = dlncg_dlnvs(1:kmaxRc,i)*thkm(1)/2.0 + !else if (i == mmax) then + ! dlncg_dlnvs(1:kmaxRc,i) = dlncg_dlnvs(1:kmaxRc,i)*thkm(mmax-1)/2.0 + !else + ! dlncg_dlnvs(1:kmaxRc,i) = dlncg_dlnvs(1:kmaxRc,i)*(thkm(i-1)+thkm(i))/2.0 + !endif + vpm(i) = vpz(i) - 0.5*dlnVp*vpz(i) call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,& @@ -110,6 +121,16 @@ subroutine depthkernel(nx,ny,nz,vel,pvRc,sen_vsRc,sen_vpRc,sen_rhoRc, & ! cga = 0.5*(cg1(nn)+cg2(nn)) dlncg_dlnvp(nn,i) = (cg2(nn)-cg1(nn))/(dlnVp*vpz(i)) enddo + + + !if (i == 1) then + ! dlncg_dlnvp(1:kmaxRc,i) = dlncg_dlnvp(1:kmaxRc,i)*thkm(1)/2.0 + !else if (i == mmax) then + ! dlncg_dlnvp(1:kmaxRc,i) = dlncg_dlnvp(1:kmaxRc,i)*thkm(mmax-1)/2.0 + !else + ! dlncg_dlnvp(1:kmaxRc,i) = dlncg_dlnvp(1:kmaxRc,i)*(thkm(i-1)+thkm(i))/2.0 + !endif + rhom(i) = rhoz(i) - 0.5*dlnrho*rhoz(i) call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,& rmax,rdep,rvp,rvs,rrho,rthk) @@ -127,6 +148,15 @@ subroutine depthkernel(nx,ny,nz,vel,pvRc,sen_vsRc,sen_vpRc,sen_rhoRc, & ! cga = 0.5*(cg1(nn)+cg2(nn)) dlncg_dlnrho(nn,i) = (cg2(nn)-cg1(nn))/(dlnrho*rhoz(i)) enddo + + ! if (i == 1) then + ! dlncg_dlnrho(1:kmaxRc,i) = dlncg_dlnrho(1:kmaxRc,i)*thkm(1)/2.0 + ! else if (i == mmax) then + ! dlncg_dlnrho(1:kmaxRc,i) = dlncg_dlnrho(1:kmaxRc,i)*thkm(mmax-1)/2.0 + ! else + ! dlncg_dlnrho(1:kmaxRc,i) = dlncg_dlnrho(1:kmaxRc,i)*(thkm(i-1)+thkm(i))/2.0 + ! endif + enddo sen_vsRc((jj-1)*nx+ii,1:kmaxRc,1:mmax)=dlncg_dlnvs(1:kmaxRc,1:mmax) sen_vpRc((jj-1)*nx+ii,1:kmaxRc,1:mmax)=dlncg_dlnvp(1:kmaxRc,1:mmax) @@ -1088,7 +1118,7 @@ subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, & iwave=1 igr=0 call depthkernel(nx,ny,nz,vels,pvLc,sen_vsLc,sen_vpLc, & - sen_rhoLc,iwave,igr,kmax,tLc,depz,minthk) + sen_rhoLc,iwave,igr,kmaxLc,tLc,depz,minthk) endif if(kmaxLg.gt.0) then diff --git a/src/DSurfTomo b/src/DSurfTomo index 59c37c3..dd8cd53 100755 Binary files a/src/DSurfTomo and b/src/DSurfTomo differ diff --git a/src/main.f90 b/src/main.f90 index 391f544..8ab9d53 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -181,9 +181,9 @@ program SurfTomo if (checkstat > 0) stop 'error allocating RP' read(10,*)(tRc(i),i=1,kmaxRc) write(*,*)'Rayleigh wave phase velocity used,periods:(s)' - write(*,'(50f6.2)')(tRc(i),i=1,kmaxRc) + write(*,'(50f7.1)')(tRc(i),i=1,kmaxRc) write(66,*)'Rayleigh wave phase velocity used,periods:(s)' - write(66,'(50f6.2)')(tRc(i),i=1,kmaxRc) + write(66,'(50f7.1)')(tRc(i),i=1,kmaxRc) endif read(10,*)kmaxRg if(kmaxRg.gt.0)then @@ -191,9 +191,9 @@ program SurfTomo if (checkstat > 0) stop 'error allocating RP' read(10,*)(tRg(i),i=1,kmaxRg) write(*,*)'Rayleigh wave group velocity used,periods:(s)' - write(*,'(50f6.2)')(tRg(i),i=1,kmaxRg) + write(*,'(50f7.1)')(tRg(i),i=1,kmaxRg) write(66,*)'Rayleigh wave group velocity used,periods:(s)' - write(66,'(50f6.2)')(tRg(i),i=1,kmaxRg) + write(66,'(50f7.1)')(tRg(i),i=1,kmaxRg) endif read(10,*)kmaxLc if(kmaxLc.gt.0)then @@ -201,9 +201,9 @@ program SurfTomo if (checkstat > 0) stop 'error allocating RP' read(10,*)(tLc(i),i=1,kmaxLc) write(*,*)'Love wave phase velocity used,periods:(s)' - write(*,'(50f6.2)')(tLc(i),i=1,kmaxLc) + write(*,'(50f7.1)')(tLc(i),i=1,kmaxLc) write(66,*)'Love wave phase velocity used,periods:(s)' - write(66,'(50f6.2)')(tLc(i),i=1,kmaxLc) + write(66,'(50f7.1)')(tLc(i),i=1,kmaxLc) endif read(10,*)kmaxLg if(kmaxLg.gt.0)then @@ -211,9 +211,9 @@ program SurfTomo if (checkstat > 0) stop 'error allocating RP' read(10,*)(tLg(i),i=1,kmaxLg) write(*,*)'Love wave group velocity used,periods:(s)' - write(*,'(50f6.2)')(tLg(i),i=1,kmaxLg) + write(*,'(50f7.1)')(tLg(i),i=1,kmaxLg) write(66,*)'Love wave group velocity used,periods:(s)' - write(66,'(50f6.2)')(tLg(i),i=1,kmaxLg) + write(66,'(50f7.1)')(tLg(i),i=1,kmaxLg) endif read(10,*)ifsyn read(10,*)noiselevel @@ -333,7 +333,7 @@ program SurfTomo enddo close(10) write(*,*) 'grid points in depth direction:(km)' - write(*,'(50f6.2)') depz + write(*,'(50f7.1)') depz @@ -411,7 +411,7 @@ program SurfTomo ! if(abs(cbst(i)) > threshold) then ! datweight(i) = exp(-(abs(cbst(i))-threshold)) ! endif - datweight(i) = 1.0/(1+0.05*exp(cbst(i)**2*threshold0)) + datweight(i) = 0.01+1.0/(1+0.05*exp(cbst(i)**2*threshold0)) cbst(i) = cbst(i)*datweight(i) enddo diff --git a/src/surfdisp96.f b/src/surfdisp96.f index 1e61103..78cc0b9 100644 --- a/src/surfdisp96.f +++ b/src/surfdisp96.f @@ -56,7 +56,7 @@ c----- integer NL, NL2, NLAY parameter(NL=200,NLAY=200,NL2=NL+NL) integer NP - parameter (NP=60) + parameter (NP=80) c----- c LIN - unit for FORTRAN read from terminal @@ -500,7 +500,7 @@ c 2 - Rayleigh Wave c iflag I*4 0 - Initialize c 1 - Make model for Love or Rayleigh Wave c----- - parameter(NL=200,NP=60) + parameter(NL=200,NP=80) real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL) integer mmax,llw c common/modl/ d,a,b,rho,rtp,dtp,btp @@ -563,7 +563,7 @@ c nev = 0 force interval halving c nev = 1 permit neville iteration if conditions are proper c nev = 2 neville iteration is being used c----- - parameter (NL=200,NP=60) + parameter (NL=200,NP=80) implicit double precision (a-h,o-z) real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL) dimension x(20),y(20) @@ -704,7 +704,7 @@ c function dltar1(wvno,omega,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) c find SH dispersion values. c - parameter (NL=200,NP=60) + parameter (NL=200,NP=80) implicit double precision (a-h,o-z) real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL) integer llw,mmax @@ -768,7 +768,7 @@ c c & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) c find P-SV dispersion values. c - parameter (NL=200,NP=60) + parameter (NL=200,NP=80) implicit double precision (a-h,o-z) dimension e(5),ee(5),ca(5,5) real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL)