diff --git a/bin/DSurfTomo b/bin/DSurfTomo new file mode 100755 index 0000000..74536d5 Binary files /dev/null and b/bin/DSurfTomo differ diff --git a/example/.gmtcommands4 b/example/.gmtcommands4 new file mode 100644 index 0000000..e3d01fa --- /dev/null +++ b/example/.gmtcommands4 @@ -0,0 +1,10 @@ +# GMT common arguments shelf +-Ba0.2f0.1:S velocity (km/s): +-JM2i +-R121.35/121.605/24.975/25.2 +-x2.3i +-Y5i +-y-2.7i +-jM2i +EOF + diff --git a/example/.gmtdefaults4 b/example/.gmtdefaults4 new file mode 100644 index 0000000..dd2074c --- /dev/null +++ b/example/.gmtdefaults4 @@ -0,0 +1,107 @@ +# +# GMT-SYSTEM 4.5.9 [64-bit] Defaults file +# +#-------- Plot Media Parameters ------------- +PAGE_COLOR = white +PAGE_ORIENTATION = landscape +PAPER_MEDIA = a4 +#-------- Basemap Annotation Parameters ------ +ANNOT_MIN_ANGLE = 20 +ANNOT_MIN_SPACING = 0 +ANNOT_FONT_PRIMARY = Helvetica +ANNOT_FONT_SIZE_PRIMARY = 6p +ANNOT_OFFSET_PRIMARY = 0.2c +ANNOT_FONT_SECONDARY = Helvetica +ANNOT_FONT_SIZE_SECONDARY = 16p +ANNOT_OFFSET_SECONDARY = 0.2c +DEGREE_SYMBOL = ring +HEADER_FONT = Helvetica +HEADER_FONT_SIZE = 36p +HEADER_OFFSET = 0.5c +LABEL_FONT = Helvetica +LABEL_FONT_SIZE = 5p +LABEL_OFFSET = 0.1c +OBLIQUE_ANNOTATION = 1 +PLOT_CLOCK_FORMAT = hh:mm:ss +PLOT_DATE_FORMAT = yyyy-mm-dd +PLOT_DEGREE_FORMAT = ddd:mm:ss +Y_AXIS_TYPE = hor_text +#-------- Basemap Layout Parameters --------- +BASEMAP_AXES = WESN +BASEMAP_FRAME_RGB = black +BASEMAP_TYPE = fancy +FRAME_PEN = 1.25p +FRAME_WIDTH = 0.1c +GRID_CROSS_SIZE_PRIMARY = 0c +GRID_PEN_PRIMARY = 0.25p +GRID_CROSS_SIZE_SECONDARY = 0c +GRID_PEN_SECONDARY = 0.5p +MAP_SCALE_HEIGHT = 0.2c +POLAR_CAP = 85/90 +TICK_LENGTH = 0.1c +TICK_PEN = 0.3p +X_AXIS_LENGTH = 25c +Y_AXIS_LENGTH = 15c +X_ORIGIN = 2.5c +Y_ORIGIN = 2.5c +UNIX_TIME = FALSE +UNIX_TIME_POS = BL/-2c/-2c +UNIX_TIME_FORMAT = %Y %b %d %H:%M:%S +#-------- Color System Parameters ----------- +COLOR_BACKGROUND = black +COLOR_FOREGROUND = white +COLOR_NAN = 128 +COLOR_IMAGE = adobe +COLOR_MODEL = rgb +HSV_MIN_SATURATION = 1 +HSV_MAX_SATURATION = 0.1 +HSV_MIN_VALUE = 0.3 +HSV_MAX_VALUE = 1 +#-------- PostScript Parameters ------------- +CHAR_ENCODING = ISOLatin1+ +DOTS_PR_INCH = 300 +GLOBAL_X_SCALE = 1 +GLOBAL_Y_SCALE = 1 +N_COPIES = 1 +PS_COLOR = rgb +PS_IMAGE_COMPRESS = lzw +PS_IMAGE_FORMAT = ascii +PS_LINE_CAP = butt +PS_LINE_JOIN = miter +PS_MITER_LIMIT = 35 +PS_VERBOSE = FALSE +TRANSPARENCY = 0 +#-------- I/O Format Parameters ------------- +D_FORMAT = %.12lg +FIELD_DELIMITER = tab +GRIDFILE_FORMAT = nf +GRIDFILE_SHORTHAND = FALSE +INPUT_CLOCK_FORMAT = hh:mm:ss +INPUT_DATE_FORMAT = yyyy-mm-dd +IO_HEADER = FALSE +N_HEADER_RECS = 1 +NAN_RECORDS = pass +OUTPUT_CLOCK_FORMAT = hh:mm:ss +OUTPUT_DATE_FORMAT = yyyy-mm-dd +OUTPUT_DEGREE_FORMAT = D +XY_TOGGLE = FALSE +#-------- Projection Parameters ------------- +ELLIPSOID = WGS-84 +MAP_SCALE_FACTOR = default +MEASURE_UNIT = cm +#-------- Calendar/Time Parameters ---------- +TIME_FORMAT_PRIMARY = full +TIME_FORMAT_SECONDARY = full +TIME_EPOCH = 2000-01-01T12:00:00 +TIME_IS_INTERVAL = OFF +TIME_INTERVAL_FRACTION = 0.5 +TIME_LANGUAGE = us +TIME_UNIT = d +TIME_WEEK_START = Sunday +Y2K_OFFSET_YEAR = 1950 +#-------- Miscellaneous Parameters ---------- +HISTORY = TRUE +INTERPOLANT = akima +LINE_STEP = 0.025c +VECTOR_SHAPE = 0 +VERBOSE = FALSE diff --git a/example/DSurfTomo.in b/example/DSurfTomo.in new file mode 100644 index 0000000..dd6984c --- /dev/null +++ b/example/DSurfTomo.in @@ -0,0 +1,21 @@ +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c INPUT PARAMETERS +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +surfdataTB.dat c: data file +18 18 9 c: nx ny nz (grid number in lat lon and depth direction) +25.2 121.35 c: goxd gozd (upper left point,[lat,lon]) +0.015 0.017 c: dvxd dvzd (grid interval in lat and lon direction) +449 c: nsrc*maxf +4.0 0.0 c: weight damp +3 c: nsublayer (numbers of sublayers for each grid interval:grid --> layer) +0.5 2.8 c: minimum velocity, maximum velocity +10 c: maxiter (iteration number) +0.2 c: sparsity fraction +26 c: kmaxRc (followed by periods) +0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 +0 c: kmaxRg +0 c: kmaxLc +0 c: kmaxLg +0 c: synthetic flag(0:real data,1:synthetic) +0.02 c: noiselevel +2.5 c: threshold diff --git a/src/CalSurfG.f90 b/src/CalSurfG.f90 new file mode 100644 index 0000000..312cfa3 --- /dev/null +++ b/src/CalSurfG.f90 @@ -0,0 +1,2838 @@ + subroutine depthkernel(nx,ny,nz,vel,pvRc,sen_vsRc,sen_vpRc,sen_rhoRc, & + iwave,igr,kmaxRc,tRc,depz,minthk) + use omp_lib + implicit none + + integer nx,ny,nz + real vel(nx,ny,nz) + real*8 sen_vpRc(ny*nx,kmaxRc,nz),sen_vsRc(ny*nx,kmaxRc,nz),sen_rhoRc(ny*nx,kmaxRc,nz) + + integer iwave,igr + real minthk + real depz(nz) + integer kmaxRc + real*8 tRc(kmaxRc) + real*8 pvRc(nx*ny,kmaxRc) + + + + real vpz(nz),vsz(nz),rhoz(nz) + real*8 dlncg_dlnvs(kmaxRc,nz),dlncg_dlnvp(kmaxRc,nz),dlncg_dlnrho(kmaxRc,nz) + integer mmax,iflsph,mode,rmax + integer ii,jj,k,i,nn,kk + integer,parameter::NL=200 + integer,parameter::NP=60 + real*8 cg1(NP),cg2(NP),cga,cgRc(NP) + real rdep(NL),rvp(NL),rvs(NL),rrho(NL),rthk(NL) + real depm(NL),vpm(NL),vsm(NL),rhom(NL),thkm(NL) + real dlnVs,dlnVp,dlnrho + + + mmax=nz + iflsph=1 + mode=1 + dlnVs=0.01 + dlnVp=0.01 + dlnrho=0.01 + + !print*,'depth kernel begin...' +!$omp parallel & +!$omp default(private) & +!$omp shared(depz,nx,ny,nz,minthk,dlnvs,dlnvp,dlnrho,kmaxRc,mmax,vel) & +!$omp shared(sen_vpRc,sen_vsRc,sen_rhoRc,tRc,pvRc,iflsph,iwave,mode,igr) +!$omp do + do jj=1,ny + do ii=1,nx + vsz(1:nz)=vel(ii,jj,1:nz) + ! some other emperical relationship maybe better, + do k=1,nz + vpz(k)=0.9409 + 2.0947*vsz(k) - 0.8206*vsz(k)**2+ & + 0.2683*vsz(k)**3 - 0.0251*vsz(k)**4 + rhoz(k)=1.6612*vpz(k) - 0.4721*vpz(k)**2 + & + 0.0671*vpz(k)**3 - 0.0043*vpz(k)**4 + & + 0.000106*vpz(k)**5 + enddo + + call refineGrid2LayerMdl(minthk,mmax,depz,vpz,vsz,rhoz,rmax,rdep,& + rvp,rvs,rrho,rthk) + call surfdisp96(rthk,rvp,rvs,rrho,rmax,iflsph,iwave,mode,igr,kmaxRc,& + tRc,cgRc) + pvRc((jj-1)*nx+ii,1:kmaxRc)=cgRc(1:kmaxRc) + !print*,cgRc(1:kmaxRc) + do kk=1,mmax-1 + depm(kk)=depz(kk) + vsm(kk) = vsz(kk) + vpm(kk) = vpz(kk) + thkm(kk) = depz(kk+1)-depz(kk) + rhom(kk) = rhoz(kk) + enddo + !!half space + depm(mmax) = depz(mmax) + vsm(mmax) = vsz(mmax) + vpm(mmax) = vpz(mmax) + rhom(mmax) = rhoz(mmax) + thkm(mmax) = 0.0 + !! calculate sensitivity kernel + do i = 1, mmax + vsm(i) = vsz(i) - 0.5*dlnVs*vsz(i) + call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,& + rmax,rdep,rvp,rvs,rrho,rthk) + call surfdisp96(rthk,rvp,rvs,rrho,rmax,iflsph,iwave,mode,& + igr,kmaxRc,tRc,cg1) + + vsm(i) = vsz(i) + 0.5*dlnVs*vsz(i) + call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,& + rmax,rdep,rvp,rvs,rrho,rthk) + call surfdisp96(rthk,rvp,rvs,rrho,rmax,iflsph,iwave,mode,& + igr,kmaxRc,tRc,cg2) + vsm(i) = vsz(i) + + do nn = 1,kmaxRc + cga = 0.5*(cg1(nn)+cg2(nn)) + dlncg_dlnvs(nn,i) = (cg2(nn)-cg1(nn))/cga/dlnVs + enddo + + + vpm(i) = vpz(i) - 0.5*dlnVp*vpz(i) + call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,& + rmax,rdep,rvp,rvs,rrho,rthk) + call surfdisp96(rthk,rvp,rvs,rrho,rmax,iflsph,iwave,mode,& + igr,kmaxRc,tRc,cg1) + + vpm(i) = vpz(i) + 0.5*dlnVp*vpz(i) + call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,& + rmax,rdep,rvp,rvs,rrho,rthk) + call surfdisp96(rthk,rvp,rvs,rrho,rmax,iflsph,iwave,mode,& + igr,kmaxRc,tRc,cg2) + vpm(i) = vpz(i) + + do nn = 1,kmaxRc + cga = 0.5*(cg1(nn)+cg2(nn)) + dlncg_dlnvp(nn,i) = (cg2(nn)-cg1(nn))/cga/dlnVp + enddo + rhom(i) = rhoz(i) - 0.5*dlnrho*rhoz(i) + call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,& + rmax,rdep,rvp,rvs,rrho,rthk) + call surfdisp96(rthk,rvp,rvs,rrho,rmax,iflsph,iwave,mode,& + igr,kmaxRc,tRc,cg1) + + rhom(i) = rhoz(i) + 0.5*dlnrho*rhoz(i) + call refineGrid2LayerMdl(minthk,mmax,depm,vpm,vsm,rhom,& + rmax,rdep,rvp,rvs,rrho,rthk) + call surfdisp96(rthk,rvp,rvs,rrho,rmax,iflsph,iwave,mode,& + igr,kmaxRc,tRc,cg2) + rhom(i) = rhoz(i) + + do nn = 1,kmaxRc + cga = 0.5*(cg1(nn)+cg2(nn)) + dlncg_dlnrho(nn,i) = (cg2(nn)-cg1(nn))/cga/dlnrho + enddo + 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) + sen_rhoRc((jj-1)*nx+ii,1:kmaxRc,1:mmax)=dlncg_dlnrho(1:kmaxRc,1:mmax) + ! print*,dlncg_dlnvp(1:kmaxRc,5) + enddo + enddo +!$omp end do +!$omp end parallel + end subroutine depthkernel + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: MODULE +! CODE: FORTRAN 90 +! This module declares variable for global use, that is, for +! USE in any subroutine or function or other module. +! Variables whose values are SAVEd can have their most +! recent values reused in any routine. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +MODULE globalp +IMPLICIT NONE +INTEGER, PARAMETER :: i10=SELECTED_REAL_KIND(6) +INTEGER :: checkstat +INTEGER, SAVE :: nvx,nvz,nnx,nnz,fom,gdx,gdz +INTEGER, SAVE :: vnl,vnr,vnt,vnb,nrnx,nrnz,sgdl,rbint +INTEGER, SAVE :: nnxr,nnzr,asgr +INTEGER, DIMENSION (:,:), ALLOCATABLE :: nsts,nstsr,srs +REAL(KIND=i10), SAVE :: gox,goz,dnx,dnz,dvx,dvz,snb,earth +REAL(KIND=i10), SAVE :: goxd,gozd,dvxd,dvzd,dnxd,dnzd +REAL(KIND=i10), SAVE :: drnx,drnz,gorx,gorz +REAL(KIND=i10), SAVE :: dnxr,dnzr,goxr,gozr +REAL(KIND=i10), DIMENSION (:,:), ALLOCATABLE, SAVE :: velv,veln,velnb +REAL(KIND=i10), DIMENSION (:,:), ALLOCATABLE, SAVE :: ttn,ttnr +!REAL(KIND=i10), DIMENSION (:), ALLOCATABLE, SAVE :: rcx,rcz +REAL(KIND=i10), PARAMETER :: pi=3.1415926535898 +!!!-------------------------------------------------------------- +!! modified by Hongjian Fang @ USTC +! real,dimension(:),allocatable,save::rw +! integer,dimension(:),allocatable,save::iw,col +! real,dimension(:,:,:),allocatable::vpf,vsf +! real,dimension(:),allocatable,save::obst,cbst,wt,dtres +!! integer,dimension(:),allocatable,save::cbst_stat +! real,dimension(:,:,:),allocatable,save::sen_vs,sen_vp,sen_rho +!!! real,dimension(:,:,:),allocatable,save::sen_vsRc,sen_vpRc,sen_rhoRc +!!! real,dimension(:,:,:),allocatable,save::sen_vsRg,sen_vpRg,sen_rhoRg +!!! real,dimension(:,:,:),allocatable,save::sen_vsLc,sen_vpLc,sen_rhoLc +!!! real,dimension(:,:,:),allocatable,save::sen_vsLg,sen_vpLg,sen_rhoLg +!!! integer,save:: count1,count2 +! integer*8,save:: nar +! integer,save:: iter,maxiter +!!!-------------------------------------------------------------- +! +! nvx,nvz = B-spline vertex values +! dvx,dvz = B-spline vertex separation +! velv(i,j) = velocity values at control points +! nnx,nnz = Number of nodes of grid in x and z +! nnxr,nnzr = Number of nodes of refined grid in x and z +! gox,goz = Origin of grid (theta,phi) +! goxr, gozr = Origin of refined grid (theta,phi) +! dnx,dnz = Node separation of grid in x and z +! dnxr,dnzr = Node separation of refined grid in x and z +! veln(i,j) = velocity values on a refined grid of nodes +! velnb(i,j) = Backup of veln required for source grid refinement +! ttn(i,j) = traveltime field on the refined grid of nodes +! ttnr(i,j) = ttn for refined grid +! nsts(i,j) = node status (-1=far,0=alive,>0=close) +! nstsr(i,j) = nsts for refined grid +! checkstat = check status of memory allocation +! fom = use first-order(0) or mixed-order(1) scheme +! snb = Maximum size of narrow band as fraction of nnx*nnz +! nrc = number of receivers +! rcx(i),rcz(i) = (x,z) coordinates of receivers +! earth = radius of Earth (in km) +! goxd,gozd = gox,goz in degrees +! dvxd,dvzd = dvx,dvz in degrees +! dnzd,dnzd = dnx,dnz in degrees +! gdx,gdz = grid dicing in x and z +! vnl,vnr,vnb,vnt = Bounds of refined grid +! nrnx,nrnz = Number of nodes in x and z for refined grid +! gorx,gorz = Grid origin of refined grid +! sgdl = Source grid dicing level +! rbint = Ray-boundary intersection (0=no, 1=yes). +! asgr = Apply source grid refinement (0=no,1=yes) +! srs = Source-receiver status (0=no path, 1=path exists) +! +END MODULE globalp + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: MODULE +! CODE: FORTRAN 90 +! This module contains all the subroutines used to calculate +! the first-arrival traveltime field through the grid. +! Subroutines are: +! (1) travel +! (2) fouds1 +! (3) fouds2 +! (4) addtree +! (5) downtree +! (6) updtree +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +MODULE traveltime +USE globalp +IMPLICIT NONE +INTEGER ntr +TYPE backpointer + INTEGER(KIND=2) :: px,pz +END TYPE backpointer +TYPE(backpointer), DIMENSION (:), ALLOCATABLE :: btg +! +! btg = backpointer to relate grid nodes to binary tree entries +! px = grid-point in x +! pz = grid-point in z +! ntr = number of entries in binary tree +! + +CONTAINS + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: SUBROUTINE +! CODE: FORTRAN 90 +! This subroutine is passed the location of a source, and from +! this point the first-arrival traveltime field through the +! velocity grid is determined. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +SUBROUTINE travel(scx,scz,urg) +IMPLICIT NONE +INTEGER :: isx,isz,sw,i,j,ix,iz,urg,swrg +REAL(KIND=i10) :: scx,scz,vsrc,dsx,dsz,ds +REAL(KIND=i10), DIMENSION (2,2) :: vss +! isx,isz = grid cell indices (i,j,k) which contains source +! scx,scz = (r,x,y) location of source +! sw = a switch (0=off,1=on) +! ix,iz = j,k position of "close" point with minimum traveltime +! maxbt = maximum size of narrow band binary tree +! rd2,rd3 = substitution variables +! vsrc = velocity at source +! vss = velocity at nodes surrounding source +! dsx, dsz = distance from source to cell boundary in x and z +! ds = distance from source to nearby node +! urg = use refined grid (0=no,1=yes,2=previously used) +! swrg = switch to end refined source grid computation +! +! The first step is to find out where the source resides +! in the grid of nodes. The cell in which it resides is +! identified by the "north-west" node of the cell. If the +! source lies on the edge or corner (a node) of the cell, then +! this scheme still applies. +! +isx=INT((scx-gox)/dnx)+1 +isz=INT((scz-goz)/dnz)+1 +sw=0 +IF(isx.lt.1.or.isx.gt.nnx)sw=1 +IF(isz.lt.1.or.isz.gt.nnz)sw=1 +IF(sw.eq.1)then + isx=90.0-isx*180.0/pi + isz=isz*180.0/pi + WRITE(6,*)"Source lies outside bounds of model (lat,long)= ",isx,isz + WRITE(6,*)"TERMINATING PROGRAM!!!" + STOP +ENDIF +IF(isx.eq.nnx)isx=isx-1 +IF(isz.eq.nnz)isz=isz-1 +! +! Set all values of nsts to -1 if beginning from a source +! point. +! +IF(urg.NE.2)nsts=-1 +! +! set initial size of binary tree to zero +! +ntr=0 +IF(urg.EQ.2)THEN +! +! In this case, source grid refinement has been applied, so +! the initial narrow band will come from resampling the +! refined grid. +! + DO i=1,nnx + DO j=1,nnz + IF(nsts(j,i).GT.0)THEN + CALL addtree(j,i) + ENDIF + ENDDO + ENDDO +ELSE +! +! In general, the source point need not lie on a grid point. +! Bi-linear interpolation is used to find velocity at the +! source point. +! + nsts=-1 + DO i=1,2 + DO j=1,2 + vss(i,j)=veln(isz-1+j,isx-1+i) + ENDDO + ENDDO + dsx=(scx-gox)-(isx-1)*dnx + dsz=(scz-goz)-(isz-1)*dnz + CALL bilinear(vss,dsx,dsz,vsrc) +! +! Now find the traveltime at the four surrounding grid points. This +! is calculated approximately by assuming the traveltime from the +! source point to each node is equal to the the distance between +! the two points divided by the average velocity of the points +! + DO i=1,2 + DO j=1,2 + ds=SQRT((dsx-(i-1)*dnx)**2+(dsz-(j-1)*dnz)**2) + ttn(isz-1+j,isx-1+i)=2.0*ds/(vss(i,j)+vsrc) + CALL addtree(isz-1+j,isx-1+i) + ENDDO + ENDDO +ENDIF +! +! Now calculate the first-arrival traveltimes at the +! remaining grid points. This is done via a loop which +! repeats the procedure of finding the first-arrival +! of all "close" points, adding it to the set of "alive" +! points and updating the points surrounding the new "alive" +! point. The process ceases when the binary tree is empty, +! in which case all grid points are "alive". +! +DO WHILE(ntr.gt.0) +! +! First, check whether source grid refinement is +! being applied; if so, then there is a special +! exit condition. +! +IF(urg.EQ.1)THEN + ix=btg(1)%px + iz=btg(1)%pz + swrg=0 + IF(ix.EQ.1)THEN + IF(vnl.NE.1)swrg=1 + ENDIF + IF(ix.EQ.nnx)THEN + IF(vnr.NE.nnx)swrg=1 + ENDIF + IF(iz.EQ.1)THEN + IF(vnt.NE.1)swrg=1 + ENDIF + IF(iz.EQ.nnz)THEN + IF(vnb.NE.nnz)swrg=1 + ENDIF + IF(swrg.EQ.1)THEN + nsts(iz,ix)=0 + EXIT + ENDIF +ENDIF +! +! Set the "close" point with minimum traveltime +! to "alive" +! + ix=btg(1)%px + iz=btg(1)%pz + nsts(iz,ix)=0 +! +! Update the binary tree by removing the root and +! sweeping down the tree. +! + CALL downtree +! +! Now update or find values of up to four grid points +! that surround the new "alive" point. +! +! Test points that vary in x +! + DO i=ix-1,ix+1,2 + IF(i.ge.1.and.i.le.nnx)THEN + IF(nsts(iz,i).eq.-1)THEN +! +! This option occurs when a far point is added to the list +! of "close" points +! + IF(fom.eq.0)THEN + CALL fouds1(iz,i) + ELSE + CALL fouds2(iz,i) + ENDIF + CALL addtree(iz,i) + ELSE IF(nsts(iz,i).gt.0)THEN +! +! This happens when a "close" point is updated +! + IF(fom.eq.0)THEN + CALL fouds1(iz,i) + ELSE + CALL fouds2(iz,i) + ENDIF + CALL updtree(iz,i) + ENDIF + ENDIF + ENDDO +! +! Test points that vary in z +! + DO i=iz-1,iz+1,2 + IF(i.ge.1.and.i.le.nnz)THEN + IF(nsts(i,ix).eq.-1)THEN +! +! This option occurs when a far point is added to the list +! of "close" points +! + IF(fom.eq.0)THEN + CALL fouds1(i,ix) + ELSE + CALL fouds2(i,ix) + ENDIF + CALL addtree(i,ix) + ELSE IF(nsts(i,ix).gt.0)THEN +! +! This happens when a "close" point is updated +! + IF(fom.eq.0)THEN + CALL fouds1(i,ix) + ELSE + CALL fouds2(i,ix) + ENDIF + CALL updtree(i,ix) + ENDIF + ENDIF + ENDDO +ENDDO +END SUBROUTINE travel + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: SUBROUTINE +! CODE: FORTRAN 90 +! This subroutine calculates a trial first-arrival traveltime +! at a given node from surrounding nodes using the +! First-Order Upwind Difference Scheme (FOUDS) of +! Sethian and Popovici (1999). +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +SUBROUTINE fouds1(iz,ix) +IMPLICIT NONE +INTEGER :: j,k,ix,iz,tsw1,swsol +REAL(KIND=i10) :: trav,travm,slown,tdsh,tref +REAL(KIND=i10) :: a,b,c,u,v,em,ri,risti +REAL(KIND=i10) :: rd1 +! +! ix = NS position of node coordinate for determination +! iz = EW vertical position of node coordinate for determination +! trav = traveltime calculated for trial node +! travm = minimum traveltime calculated for trial node +! slown = slowness at (iz,ix) +! tsw1 = traveltime switch (0=first time,1=previously) +! a,b,c,u,v,em = Convenience variables for solving quadratic +! tdsh = local traveltime from neighbouring node +! tref = reference traveltime at neighbouring node +! ri = Radial distance +! risti = ri*sin(theta) at point (iz,ix) +! rd1 = dummy variable +! swsol = switch for solution (0=no solution, 1=solution) +! +! Inspect each of the four quadrants for the minimum time +! solution. +! +tsw1=0 +slown=1.0/veln(iz,ix) +ri=earth +risti=ri*sin(gox+(ix-1)*dnx) +DO j=ix-1,ix+1,2 + DO k=iz-1,iz+1,2 + IF(j.GE.1.AND.j.LE.nnx)THEN + IF(k.GE.1.AND.k.LE.nnz)THEN +! +! There are seven solution options in +! each quadrant. +! + swsol=0 + IF(nsts(iz,j).EQ.0)THEN + swsol=1 + IF(nsts(k,ix).EQ.0)THEN + u=ri*dnx + v=risti*dnz + em=ttn(k,ix)-ttn(iz,j) + a=u**2+v**2 + b=-2.0*u**2*em + c=u**2*(em**2-v**2*slown**2) + tref=ttn(iz,j) + ELSE + a=1.0 + b=0.0 + c=-slown**2*ri**2*dnx**2 + tref=ttn(iz,j) + ENDIF + ELSE IF(nsts(k,ix).EQ.0)THEN + swsol=1 + a=1.0 + b=0.0 + c=-(slown*risti*dnz)**2 + tref=ttn(k,ix) + ENDIF +! +! Now find the solution of the quadratic equation +! + IF(swsol.EQ.1)THEN + rd1=b**2-4.0*a*c + IF(rd1.LT.0.0)rd1=0.0 + tdsh=(-b+sqrt(rd1))/(2.0*a) + trav=tref+tdsh + IF(tsw1.EQ.1)THEN + travm=MIN(trav,travm) + ELSE + travm=trav + tsw1=1 + ENDIF + ENDIF + ENDIF + ENDIF + ENDDO +ENDDO +ttn(iz,ix)=travm +END SUBROUTINE fouds1 + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: SUBROUTINE +! CODE: FORTRAN 90 +! This subroutine calculates a trial first-arrival traveltime +! at a given node from surrounding nodes using the +! Mixed-Order (2nd) Upwind Difference Scheme (FOUDS) of +! Popovici and Sethian (2002). +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +SUBROUTINE fouds2(iz,ix) +IMPLICIT NONE +INTEGER :: j,k,j2,k2,ix,iz,tsw1 +INTEGER :: swj,swk,swsol +REAL(KIND=i10) :: trav,travm,slown,tdsh,tref,tdiv +REAL(KIND=i10) :: a,b,c,u,v,em,ri,risti,rd1 +! +! ix = NS position of node coordinate for determination +! iz = EW vertical position of node coordinate for determination +! trav = traveltime calculated for trial node +! travm = minimum traveltime calculated for trial node +! slown = slowness at (iz,ix) +! tsw1 = traveltime switch (0=first time,1=previously) +! a,b,c,u,v,em = Convenience variables for solving quadratic +! tdsh = local traveltime from neighbouring node +! tref = reference traveltime at neighbouring node +! ri = Radial distance +! risti = ri*sin(theta) at point (iz,ix) +! swj,swk = switches for second order operators +! tdiv = term to divide tref by depending on operator order +! swsol = switch for solution (0=no solution, 1=solution) +! +! Inspect each of the four quadrants for the minimum time +! solution. +! +tsw1=0 +slown=1.0/veln(iz,ix) +ri=earth +risti=ri*sin(gox+(ix-1)*dnx) +DO j=ix-1,ix+1,2 + IF(j.GE.1.AND.j.LE.nnx)THEN + swj=-1 + IF(j.eq.ix-1)THEN + j2=j-1 + IF(j2.GE.1)THEN + IF(nsts(iz,j2).EQ.0)swj=0 + ENDIF + ELSE + j2=j+1 + IF(j2.LE.nnx)THEN + IF(nsts(iz,j2).EQ.0)swj=0 + ENDIF + ENDIF + IF(nsts(iz,j).EQ.0.AND.swj.EQ.0)THEN + swj=-1 + IF(ttn(iz,j).GT.ttn(iz,j2))THEN + swj=0 + ENDIF + ELSE + swj=-1 + ENDIF + DO k=iz-1,iz+1,2 + IF(k.GE.1.AND.k.LE.nnz)THEN + swk=-1 + IF(k.eq.iz-1)THEN + k2=k-1 + IF(k2.GE.1)THEN + IF(nsts(k2,ix).EQ.0)swk=0 + ENDIF + ELSE + k2=k+1 + IF(k2.LE.nnz)THEN + IF(nsts(k2,ix).EQ.0)swk=0 + ENDIF + ENDIF + IF(nsts(k,ix).EQ.0.AND.swk.EQ.0)THEN + swk=-1 + IF(ttn(k,ix).GT.ttn(k2,ix))THEN + swk=0 + ENDIF + ELSE + swk=-1 + ENDIF +! +! There are 8 solution options in +! each quadrant. +! + swsol=0 + IF(swj.EQ.0)THEN + swsol=1 + IF(swk.EQ.0)THEN + u=2.0*ri*dnx + v=2.0*risti*dnz + em=4.0*ttn(iz,j)-ttn(iz,j2)-4.0*ttn(k,ix) + em=em+ttn(k2,ix) + a=v**2+u**2 + b=2.0*em*u**2 + c=u**2*(em**2-slown**2*v**2) + tref=4.0*ttn(iz,j)-ttn(iz,j2) + tdiv=3.0 + ELSE IF(nsts(k,ix).EQ.0)THEN + u=risti*dnz + v=2.0*ri*dnx + em=3.0*ttn(k,ix)-4.0*ttn(iz,j)+ttn(iz,j2) + a=v**2+9.0*u**2 + b=6.0*em*u**2 + c=u**2*(em**2-slown**2*v**2) + tref=ttn(k,ix) + tdiv=1.0 + ELSE + u=2.0*ri*dnx + a=1.0 + b=0.0 + c=-u**2*slown**2 + tref=4.0*ttn(iz,j)-ttn(iz,j2) + tdiv=3.0 + ENDIF + ELSE IF(nsts(iz,j).EQ.0)THEN + swsol=1 + IF(swk.EQ.0)THEN + u=ri*dnx + v=2.0*risti*dnz + em=3.0*ttn(iz,j)-4.0*ttn(k,ix)+ttn(k2,ix) + a=v**2+9.0*u**2 + b=6.0*em*u**2 + c=u**2*(em**2-v**2*slown**2) + tref=ttn(iz,j) + tdiv=1.0 + ELSE IF(nsts(k,ix).EQ.0)THEN + u=ri*dnx + v=risti*dnz + em=ttn(k,ix)-ttn(iz,j) + a=u**2+v**2 + b=-2.0*u**2*em + c=u**2*(em**2-v**2*slown**2) + tref=ttn(iz,j) + tdiv=1.0 + ELSE + a=1.0 + b=0.0 + c=-slown**2*ri**2*dnx**2 + tref=ttn(iz,j) + tdiv=1.0 + ENDIF + ELSE + IF(swk.EQ.0)THEN + swsol=1 + u=2.0*risti*dnz + a=1.0 + b=0.0 + c=-u**2*slown**2 + tref=4.0*ttn(k,ix)-ttn(k2,ix) + tdiv=3.0 + ELSE IF(nsts(k,ix).EQ.0)THEN + swsol=1 + a=1.0 + b=0.0 + c=-slown**2*risti**2*dnz**2 + tref=ttn(k,ix) + tdiv=1.0 + ENDIF + ENDIF +! +! Now find the solution of the quadratic equation +! + IF(swsol.EQ.1)THEN + rd1=b**2-4.0*a*c + IF(rd1.LT.0.0)rd1=0.0 + tdsh=(-b+sqrt(rd1))/(2.0*a) + trav=(tref+tdsh)/tdiv + IF(tsw1.EQ.1)THEN + travm=MIN(trav,travm) + ELSE + travm=trav + tsw1=1 + ENDIF + ENDIF + ENDIF + ENDDO + ENDIF +ENDDO +ttn(iz,ix)=travm +END SUBROUTINE fouds2 + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: SUBROUTINE +! CODE: FORTRAN 90 +! This subroutine adds a value to the binary tree by +! placing a value at the bottom and pushing it up +! to its correct position. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +SUBROUTINE addtree(iz,ix) +IMPLICIT NONE +INTEGER :: ix,iz,tpp,tpc +TYPE(backpointer) :: exch +! +! ix,iz = grid position of new addition to tree +! tpp = tree position of parent +! tpc = tree position of child +! exch = dummy to exchange btg values +! +! First, increase the size of the tree by one. +! +ntr=ntr+1 +! +! Put new value at base of tree +! +nsts(iz,ix)=ntr +btg(ntr)%px=ix +btg(ntr)%pz=iz +! +! Now filter the new value up to its correct position +! +tpc=ntr +tpp=tpc/2 +DO WHILE(tpp.gt.0) + IF(ttn(iz,ix).lt.ttn(btg(tpp)%pz,btg(tpp)%px))THEN + nsts(iz,ix)=tpp + nsts(btg(tpp)%pz,btg(tpp)%px)=tpc + exch=btg(tpc) + btg(tpc)=btg(tpp) + btg(tpp)=exch + tpc=tpp + tpp=tpc/2 + ELSE + tpp=0 + ENDIF +ENDDO +END SUBROUTINE addtree + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: SUBROUTINE +! CODE: FORTRAN 90 +! This subroutine updates the binary tree after the root +! value has been used. The root is replaced by the value +! at the bottom of the tree, which is then filtered down +! to its correct position. This ensures that the tree remains +! balanced. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +SUBROUTINE downtree +IMPLICIT NONE +INTEGER :: tpp,tpc +REAL(KIND=i10) :: rd1,rd2 +TYPE(backpointer) :: exch +! +! tpp = tree position of parent +! tpc = tree position of child +! exch = dummy to exchange btg values +! rd1,rd2 = substitution variables +! +! Replace root of tree with its last value +! +IF(ntr.EQ.1)THEN + ntr=ntr-1 + RETURN +ENDIF +nsts(btg(ntr)%pz,btg(ntr)%px)=1 +btg(1)=btg(ntr) +! +! Reduce size of tree by one +! +ntr=ntr-1 +! +! Now filter new root down to its correct position +! +tpp=1 +tpc=2*tpp +DO WHILE(tpc.lt.ntr) +! +! Check which of the two children is smallest - use the smallest +! + rd1=ttn(btg(tpc)%pz,btg(tpc)%px) + rd2=ttn(btg(tpc+1)%pz,btg(tpc+1)%px) + IF(rd1.gt.rd2)THEN + tpc=tpc+1 + ENDIF +! +! Check whether the child is smaller than the parent; if so, then swap, +! if not, then we are done +! + rd1=ttn(btg(tpc)%pz,btg(tpc)%px) + rd2=ttn(btg(tpp)%pz,btg(tpp)%px) + IF(rd1.lt.rd2)THEN + nsts(btg(tpp)%pz,btg(tpp)%px)=tpc + nsts(btg(tpc)%pz,btg(tpc)%px)=tpp + exch=btg(tpc) + btg(tpc)=btg(tpp) + btg(tpp)=exch + tpp=tpc + tpc=2*tpp + ELSE + tpc=ntr+1 + ENDIF +ENDDO +! +! If ntr is an even number, then we still have one more test to do +! +IF(tpc.eq.ntr)THEN + rd1=ttn(btg(tpc)%pz,btg(tpc)%px) + rd2=ttn(btg(tpp)%pz,btg(tpp)%px) + IF(rd1.lt.rd2)THEN + nsts(btg(tpp)%pz,btg(tpp)%px)=tpc + nsts(btg(tpc)%pz,btg(tpc)%px)=tpp + exch=btg(tpc) + btg(tpc)=btg(tpp) + btg(tpp)=exch + ENDIF +ENDIF +END SUBROUTINE downtree + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: SUBROUTINE +! CODE: FORTRAN 90 +! This subroutine updates a value on the binary tree. The FMM +! should only produce updated values that are less than their +! prior values. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +SUBROUTINE updtree(iz,ix) +IMPLICIT NONE +INTEGER :: ix,iz,tpp,tpc +TYPE(backpointer) :: exch +! +! ix,iz = grid position of new addition to tree +! tpp = tree position of parent +! tpc = tree position of child +! exch = dummy to exchange btg values +! +! Filter the updated value to its correct position +! +tpc=nsts(iz,ix) +tpp=tpc/2 +DO WHILE(tpp.gt.0) + IF(ttn(iz,ix).lt.ttn(btg(tpp)%pz,btg(tpp)%px))THEN + nsts(iz,ix)=tpp + nsts(btg(tpp)%pz,btg(tpp)%px)=tpc + exch=btg(tpc) + btg(tpc)=btg(tpp) + btg(tpp)=exch + tpc=tpp + tpp=tpc/2 + ELSE + tpp=0 + ENDIF +ENDDO +END SUBROUTINE updtree + +END MODULE traveltime + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! MAIN PROGRAM +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: PROGRAM +! CODE: FORTRAN 90 +! This program is designed to implement the Fast Marching +! Method (FMM) for calculating first-arrival traveltimes +! through a 2-D continuous velocity medium in spherical shell +! coordinates (x=theta or latitude, z=phi or longitude). +! It is written in Fortran 90, although it is probably more +! accurately described as Fortran 77 with some of the Fortran 90 +! extensions. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!PROGRAM tomo_surf +subroutine CalSurfG(nx,ny,nz,nparpi,vels,iw,rw,col,dsurf, & + goxdf,gozdf,dvxdf,dvzdf,kmaxRc,kmaxRg,kmaxLc,kmaxLg, & + tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk, & + scxf,sczf,rcxf,rczf,nrc1,nsrcsurf1,knum1,kmax,nsrcsurf,nrcf, & + nar,writepath) +USE globalp +USE traveltime +IMPLICIT NONE +!CHARACTER (LEN=30) ::grid,frechet +!CHARACTER (LEN=40) :: sources,receivers,otimes +!CHARACTER (LEN=30) :: travelt,rtravel,wrays,cdum +INTEGER :: i,j,k,l,nsrc,tnr,urg +INTEGER :: sgs,isx,isz,sw,idm1,idm2,nnxb,nnzb +INTEGER :: ogx,ogz,grdfx,grdfz,maxbt +REAL(KIND=i10) :: x,z,goxb,gozb,dnxb,dnzb +!REAL(KIND=i10), DIMENSION (:,:), ALLOCATABLE :: scxf,sczf +!REAL(KIND=i10), DIMENSION (:,:,:), ALLOCATABLE :: rcxf,rczf +! +! sources = File containing source locations +! receivers = File containing receiver locations +! grid = File containing grid of velocity vertices for +! resampling on a finer grid with cubic B-splines +! frechet = output file containing matrix of frechet derivatives +! travelt = File name for storage of traveltime field +! wttf = Write traveltimes to file? (0=no,>0=source id) +! fom = Use first-order(0) or mixed-order(1) scheme +! nsrc = number of sources +! scx,scz = source location in r,x,z +! scx,scz = source location in r,x,z +! x,z = temporary variables for source location +! fsrt = find source-receiver traveltimes? (0=no,1=yes) +! rtravel = output file for source-receiver traveltimes +! cdum = dummy character variable ! wrgf = write ray geometries to file? (<0=all,0=no,>0=source id.) +! wrays = file containing raypath geometries +! cfd = calculate Frechet derivatives? (0=no, 1=yes) +! tnr = total number of receivers +! sgs = Extent of refined source grid +! isx,isz = cell containing source +! nnxb,nnzb = Backup for nnz,nnx +! goxb,gozb = Backup for gox,goz +! dnxb,dnzb = Backup for dnx,dnz +! ogx,ogz = Location of refined grid origin +! gridfx,grdfz = Number of refined nodes per cell +! urg = use refined grid (0=no,1=yes,2=previously used) +! maxbt = maximum size of narrow band binary tree +! otimes = file containing source-receiver association information +!c----------------------------------------------------------------- +! variables defined by Hongjian Fang + integer nx,ny,nz + integer kmax,nsrcsurf,nrcf + real vels(nx,ny,nz) + real rw(*) + integer iw(*),col(*) + real dsurf(*) + real goxdf,gozdf,dvxdf,dvzdf + integer kmaxRc,kmaxRg,kmaxLc,kmaxLg + real*8 tRc(*),tRg(*),tLc(*),tLg(*) + integer wavetype(nsrcsurf,kmax) + integer periods(nsrcsurf,kmax),nrc1(nsrcsurf,kmax),nsrcsurf1(kmax) + integer knum1(kmax),igrt(nsrcsurf,kmax) + real scxf(nsrcsurf,kmax),sczf(nsrcsurf,kmax),rcxf(nrcf,nsrcsurf,kmax),rczf(nrcf,nsrcsurf,kmax) + integer nar + real minthk + integer nparpi + + + real vpz(nz),vsz(nz),rhoz(nz),depz(nz) + real*8 pvRc(nx*ny,kmaxRc),pvRg(nx*ny,kmaxRg),pvLc(nx*ny,kmaxLc),pvLg(nx*ny,kmaxLg) + real*8 sen_vsRc(nx*ny,kmaxRc,nz),sen_vpRc(nx*ny,kmaxRc,nz) + real*8 sen_rhoRc(nx*ny,kmaxRc,nz) + real*8 sen_vsRg(nx*ny,kmaxRg,nz),sen_vpRg(nx*ny,kmaxRg,nz) + real*8 sen_rhoRg(nx*ny,kmaxRg,nz) + real*8 sen_vsLc(nx*ny,kmaxLc,nz),sen_vpLc(nx*ny,kmaxLc,nz) + real*8 sen_rhoLc(nx*ny,kmaxLc,nz) + real*8 sen_vsLg(nx*ny,kmaxLg,nz),sen_vpLg(nx*ny,kmaxLg,nz) + real*8 sen_rhoLg(nx*ny,kmaxLg,nz) + real*8 sen_vs(nx*ny,kmax,nz),sen_vp(nx*ny,kmax,nz) + real*8 sen_rho(nx*ny,kmax,nz) + real coe_rho(nz-1),coe_a(nz-1) + real*8 velf(ny*nx) + integer kmax1,kmax2,kmax3,count1 + integer igr + integer iwave + integer knumi,srcnum + real,dimension(:,:),allocatable:: fdm + real row(nparpi) + real vpft(nz-1) + real cbst1 + integer ii,jj,kk,nn,istep + integer level,maxlevel,maxleveld,HorizonType,VerticalType,PorS + real,parameter::ftol=1e-4 + integer writepath +gdx=5 +gdz=5 +asgr=1 +sgdl=8 +sgs=8 +earth=6371.0 +fom=1 +snb=0.5 +goxd=goxdf +gozd=gozdf +dvxd=dvxdf +dvzd=dvzdf +nvx=nx-2 +nvz=ny-2 +ALLOCATE(velv(0:nvz+1,0:nvx+1), STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with ALLOCATE: SUBROUTINE gridder: REAL velv' +ENDIF +! +! Convert from degrees to radians +! +dvx=dvxd*pi/180.0 +dvz=dvzd*pi/180.0 +gox=(90.0-goxd)*pi/180.0 +goz=gozd*pi/180.0 +! +! Compute corresponding values for propagation grid. +! +nnx=(nvx-1)*gdx+1 +nnz=(nvz-1)*gdz+1 +dnx=dvx/gdx +dnz=dvz/gdz +dnxd=dvxd/gdx +dnzd=dvzd/gdz +ALLOCATE(veln(nnz,nnx), STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with ALLOCATE: SUBROUTINE gridder: REAL veln' +ENDIF + +! +! Call a subroutine which reads in the velocity grid +! +!CALL gridder(grid) +! +! Read in all source coordinates. +! +! +! Now work out, source by source, the first-arrival traveltime +! field plus source-receiver traveltimes +! and ray paths if required. First, allocate memory to the +! traveltime field array +! +ALLOCATE(ttn(nnz,nnx), STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with ALLOCATE: PROGRAM fmmin2d: REAL ttn' +ENDIF + rbint=0 +! +! Allocate memory for node status and binary trees +! +ALLOCATE(nsts(nnz,nnx)) +maxbt=NINT(snb*nnx*nnz) +ALLOCATE(btg(maxbt)) + +allocate(fdm(0:nvz+1,0:nvx+1)) + + if(kmaxRc.gt.0) then + iwave=2 + igr=0 + call depthkernel(nx,ny,nz,vels,pvRc,sen_vsRc,sen_vpRc, & + sen_rhoRc,iwave,igr,kmaxRc,tRc,depz,minthk) + endif + + if(kmaxRg.gt.0) then + iwave=2 + igr=1 + call depthkernel(nx,ny,nz,vels,pvRg,sen_vsRg,sen_vpRg, & + sen_rhoRg,iwave,igr,kmaxRg,tRg,depz,minthk) + endif + + if(kmaxLc.gt.0) then + iwave=1 + igr=0 + call depthkernel(nx,ny,nz,vels,pvLc,sen_vsLc,sen_vpLc, & + sen_rhoLc,iwave,igr,kmaxLc,tLc,depz,minthk) + endif + + if(kmaxLg.gt.0) then + iwave=1 + igr=1 + call depthkernel(nx,ny,nz,vels,pvLg,sen_vsLg,sen_vpLg, & + sen_rhoLg,iwave,igr,kmaxLg,tLg,depz,minthk) + endif + +nar=0 +count1=0 + +sen_vs=0 +sen_vp=0 +sen_rho=0 +kmax1=kmaxRc +kmax2=kmaxRc+kmaxRg +kmax3=kmaxRc+kmaxRg+kmaxLc +do knumi=1,kmax +do srcnum=1,nsrcsurf1(knum1(knumi)) + if(wavetype(srcnum,knum1(knumi))==2.and.igrt(srcnum,knum1(knumi))==0) then + velf(1:nx*ny)=pvRc(1:nx*ny,periods(srcnum,knum1(knumi))) + sen_vs(:,1:kmax1,:)=sen_vsRc(:,1:kmaxRc,:)!(:,nt(istep),:) + sen_vp(:,1:kmax1,:)=sen_vpRc(:,1:kmaxRc,:)!(:,nt(istep),:) + sen_rho(:,1:kmax1,:)=sen_rhoRc(:,1:kmaxRc,:)!(:,nt(istep),:) + endif + if(wavetype(srcnum,knum1(knumi))==2.and.igrt(srcnum,knum1(knumi))==1) then + velf(1:nx*ny)=pvRg(1:nx*ny,periods(srcnum,knum1(knumi))) + sen_vs(:,kmax1+1:kmax2,:)=sen_vsRg(:,1:kmaxRg,:)!(:,nt,:) + sen_vp(:,kmax1+1:kmax2,:)=sen_vpRg(:,1:kmaxRg,:)!(:,nt,:) + sen_rho(:,kmax1+1:kmax2,:)=sen_rhoRg(:,1:kmaxRg,:)!(:,nt,:) + endif + if(wavetype(srcnum,knum1(knumi))==1.and.igrt(srcnum,knum1(knumi))==0) then + velf(1:nx*ny)=pvLc(1:nx*ny,periods(srcnum,knum1(knumi))) + sen_vs(:,kmax2+1:kmax3,:)=sen_vsLc(:,1:kmaxLc,:)!(:,nt,:) + sen_vp(:,kmax2+1:kmax3,:)=sen_vpLc(:,1:kmaxLc,:)!(:,nt,:) + sen_rho(:,kmax2+1:kmax3,:)=sen_rhoLc(:,1:kmaxLc,:)!(:,nt,:) + endif + if(wavetype(srcnum,knum1(knumi))==1.and.igrt(srcnum,knum1(knumi))==1) then + velf(1:nx*ny)=pvLg(1:nx*ny,periods(srcnum,knum1(knumi))) + sen_vs(:,kmax3+1:kmax,:)=sen_vsLg(:,1:kmaxLg,:)!(:,nt,:) + sen_vp(:,kmax3+1:kmax,:)=sen_vpLg(:,1:kmaxLg,:)!(:,nt,:) + sen_rho(:,kmax3+1:kmax,:)=sen_rhoLg(:,1:kmaxLg,:)!(:,nt,:) + endif + +call gridder(velf) + x=scxf(srcnum,knum1(knumi)) + z=sczf(srcnum,knum1(knumi)) +! +! Begin by computing refined source grid if required +! + urg=0 + IF(asgr.EQ.1)THEN +! +! Back up coarse velocity grid to a holding matrix +! + ALLOCATE(velnb(nnz,nnx)) + ! MODIFIEDY BY HONGJIAN FANG @ USTC 2014/04/17 + velnb(1:nnz,1:nnx)=veln(1:nnz,1:nnx) + nnxb=nnx + nnzb=nnz + dnxb=dnx + dnzb=dnz + goxb=gox + gozb=goz +! +! Identify nearest neighbouring node to source +! + isx=INT((x-gox)/dnx)+1 + isz=INT((z-goz)/dnz)+1 + sw=0 + IF(isx.lt.1.or.isx.gt.nnx)sw=1 + IF(isz.lt.1.or.isz.gt.nnz)sw=1 + IF(sw.eq.1)then + isx=90.0-isx*180.0/pi + isz=isz*180.0/pi + WRITE(6,*)"Source lies outside bounds of model (lat,long)= ",isx,isz + WRITE(6,*)"TERMINATING PROGRAM!!!" + STOP + ENDIF + IF(isx.eq.nnx)isx=isx-1 + IF(isz.eq.nnz)isz=isz-1 +! +! Now find rectangular box that extends outward from the nearest source node +! to "sgs" nodes away. +! + vnl=isx-sgs + IF(vnl.lt.1)vnl=1 + vnr=isx+sgs + IF(vnr.gt.nnx)vnr=nnx + vnt=isz-sgs + IF(vnt.lt.1)vnt=1 + vnb=isz+sgs + IF(vnb.gt.nnz)vnb=nnz + nrnx=(vnr-vnl)*sgdl+1 + nrnz=(vnb-vnt)*sgdl+1 + drnx=dvx/REAL(gdx*sgdl) + drnz=dvz/REAL(gdz*sgdl) + gorx=gox+dnx*(vnl-1) + gorz=goz+dnz*(vnt-1) + nnx=nrnx + nnz=nrnz + dnx=drnx + dnz=drnz + gox=gorx + goz=gorz +! +! Reallocate velocity and traveltime arrays if nnx>nnxb or +! nnz 0)THEN + WRITE(6,*)'Error with DEALLOCATE: PROGRAM fmmin2d: velnb' + ENDIF +ENDIF +enddo +enddo +deallocate(fdm) +deallocate(velv,veln,ttn,nsts,btg) +END subroutine +SUBROUTINE gridder(pv) +!subroutine gridder(pv) +!subroutine gridder() +USE globalp +IMPLICIT NONE +INTEGER :: i,j,l,m,i1,j1,conx,conz,stx,stz +REAL(KIND=i10) :: u,sumi,sumj +REAL(KIND=i10), DIMENSION(:,:), ALLOCATABLE :: ui,vi +!CHARACTER (LEN=30) :: grid +! +! u = independent parameter for b-spline +! ui,vi = bspline basis functions +! conx,conz = variables for edge of B-spline grid +! stx,stz = counters for veln grid points +! sumi,sumj = summation variables for computing b-spline +! +!C--------------------------------------------------------------- +double precision pv(*) +!integer count1 +!C--------------------------------------------------------------- +! Open the grid file and read in the velocity grid. +! +!OPEN(UNIT=10,FILE=grid,STATUS='old') +!READ(10,*)nvx,nvz +!READ(10,*)goxd,gozd +!READ(10,*)dvxd,dvzd +!count1=0 +DO i=0,nvz+1 + DO j=0,nvx+1 +! count1=count1+1 +! READ(10,*)velv(i,j) +! velv(i,j)=real(pv(count1)) + velv(i,j)=real(pv(i*(nvx+2)+j+1)) + ENDDO +ENDDO +!CLOSE(10) +! +! Convert from degrees to radians +! +! +! Now dice up the grid +! +ALLOCATE(ui(gdx+1,4), STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with ALLOCATE: Subroutine gridder: REAL ui' +ENDIF +DO i=1,gdx+1 + u=gdx + u=(i-1)/u + ui(i,1)=(1.0-u)**3/6.0 + ui(i,2)=(4.0-6.0*u**2+3.0*u**3)/6.0 + ui(i,3)=(1.0+3.0*u+3.0*u**2-3.0*u**3)/6.0 + ui(i,4)=u**3/6.0 +ENDDO +ALLOCATE(vi(gdz+1,4), STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with ALLOCATE: Subroutine gridder: REAL vi' +ENDIF +DO i=1,gdz+1 + u=gdz + u=(i-1)/u + vi(i,1)=(1.0-u)**3/6.0 + vi(i,2)=(4.0-6.0*u**2+3.0*u**3)/6.0 + vi(i,3)=(1.0+3.0*u+3.0*u**2-3.0*u**3)/6.0 + vi(i,4)=u**3/6.0 +ENDDO +DO i=1,nvz-1 + conz=gdz + IF(i==nvz-1)conz=gdz+1 + DO j=1,nvx-1 + conx=gdx + IF(j==nvx-1)conx=gdx+1 + DO l=1,conz + stz=gdz*(i-1)+l + DO m=1,conx + stx=gdx*(j-1)+m + sumi=0.0 + DO i1=1,4 + sumj=0.0 + DO j1=1,4 + sumj=sumj+ui(m,j1)*velv(i-2+i1,j-2+j1) + ENDDO + sumi=sumi+vi(l,i1)*sumj + ENDDO + veln(stz,stx)=sumi + ENDDO + ENDDO + ENDDO +ENDDO +DEALLOCATE(ui,vi, STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with DEALLOCATE: SUBROUTINE gridder: REAL ui,vi' +ENDIF +END SUBROUTINE gridder + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: SUBROUTINE +! CODE: FORTRAN 90 +! This subroutine is similar to bsplreg except that it has been +! modified to deal with source grid refinement +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +SUBROUTINE bsplrefine +USE globalp +INTEGER :: i,j,k,l,i1,j1,st1,st2,nrzr,nrxr +INTEGER :: origx,origz,conx,conz,idm1,idm2 +REAL(KIND=i10) :: u,v +REAL(KIND=i10), DIMENSION (4) :: sum +REAL(KIND=i10), DIMENSION(gdx*sgdl+1,gdz*sgdl+1,4) :: ui,vi +! +! nrxr,nrzr = grid refinement level for source grid in x,z +! origx,origz = local origin of refined source grid +! +! Begin by calculating the values of the basis functions +! +nrxr=gdx*sgdl +nrzr=gdz*sgdl +DO i=1,nrzr+1 + v=nrzr + v=(i-1)/v + DO j=1,nrxr+1 + u=nrxr + u=(j-1)/u + ui(j,i,1)=(1.0-u)**3/6.0 + ui(j,i,2)=(4.0-6.0*u**2+3.0*u**3)/6.0 + ui(j,i,3)=(1.0+3.0*u+3.0*u**2-3.0*u**3)/6.0 + ui(j,i,4)=u**3/6.0 + vi(j,i,1)=(1.0-v)**3/6.0 + vi(j,i,2)=(4.0-6.0*v**2+3.0*v**3)/6.0 + vi(j,i,3)=(1.0+3.0*v+3.0*v**2-3.0*v**3)/6.0 + vi(j,i,4)=v**3/6.0 + ENDDO +ENDDO +! +! Calculate the velocity values. +! +origx=(vnl-1)*sgdl+1 +origz=(vnt-1)*sgdl+1 +DO i=1,nvz-1 + conz=nrzr + IF(i==nvz-1)conz=nrzr+1 + DO j=1,nvx-1 + conx=nrxr + IF(j==nvx-1)conx=nrxr+1 + DO k=1,conz + st1=gdz*(i-1)+(k-1)/sgdl+1 + IF(st1.LT.vnt.OR.st1.GT.vnb)CYCLE + st1=nrzr*(i-1)+k + DO l=1,conx + st2=gdx*(j-1)+(l-1)/sgdl+1 + IF(st2.LT.vnl.OR.st2.GT.vnr)CYCLE + st2=nrxr*(j-1)+l + DO i1=1,4 + sum(i1)=0.0 + DO j1=1,4 + sum(i1)=sum(i1)+ui(l,k,j1)*velv(i-2+i1,j-2+j1) + ENDDO + sum(i1)=vi(l,k,i1)*sum(i1) + ENDDO + idm1=st1-origz+1 + idm2=st2-origx+1 + IF(idm1.LT.1.OR.idm1.GT.nnz)CYCLE + IF(idm2.LT.1.OR.idm2.GT.nnx)CYCLE + veln(idm1,idm2)=sum(1)+sum(2)+sum(3)+sum(4) + ENDDO + ENDDO + ENDDO +ENDDO +END SUBROUTINE bsplrefine +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: SUBROUTINE +! CODE: FORTRAN 90 +! This subroutine calculates all receiver traveltimes for +! a given source and writes the results to file. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!SUBROUTINE srtimes(scx,scz,rcx1,rcz1,cbst1) +SUBROUTINE srtimes(scx,scz,rcx1,rcz1,cbst1) +USE globalp +IMPLICIT NONE +INTEGER :: i,k,l,irx,irz,sw,isx,isz,csid +INTEGER, PARAMETER :: noray=0,yesray=1 +INTEGER, PARAMETER :: i5=SELECTED_REAL_KIND(6) +REAL(KIND=i5) :: trr +REAL(KIND=i5), PARAMETER :: norayt=0.0 +REAL(KIND=i10) :: drx,drz,produ,scx,scz +REAL(KIND=i10) :: rcx1,rcz1,cbst1 +REAL(KIND=i10) :: sred,dpl,rd1,vels,velr +REAL(KIND=i10), DIMENSION (2,2) :: vss +!!------------------------------------------------------ +! modified by Hongjian Fang @ USTC + integer no_p,nsrc + real dist +! real cbst(*) !note that the type difference(kind=i5 vs real) +! integer cbst_stat(*) +!!------------------------------------------------------ +! +! irx,irz = Coordinates of cell containing receiver +! trr = traveltime value at receiver +! produ = dummy multiplier +! drx,drz = receiver distance from (i,j,k) grid node +! scx,scz = source coordinates +! isx,isz = source cell location +! sred = Distance from source to receiver +! dpl = Minimum path length in source neighbourhood. +! vels,velr = velocity at source and receiver +! vss = velocity at four grid points about source or receiver. +! csid = current source ID +! noray = switch to indicate no ray present +! norayt = default value given to null ray +! yesray = switch to indicate that ray is present +! +! Determine source-receiver traveltimes one at a time. +! +!0605DO i=1,nrc +!0605 IF(srs(i,csid).EQ.0)THEN +!0605! WRITE(10,*)noray,norayt +!0605 CYCLE +!0605 ENDIF +! +! The first step is to locate the receiver in the grid. +! + irx=INT((rcx1-gox)/dnx)+1 + irz=INT((rcz1-goz)/dnz)+1 + sw=0 + IF(irx.lt.1.or.irx.gt.nnx)sw=1 + IF(irz.lt.1.or.irz.gt.nnz)sw=1 + IF(sw.eq.1)then + irx=90.0-irx*180.0/pi + irz=irz*180.0/pi + WRITE(6,*)"srtimes Receiver lies outside model (lat,long)= ",irx,irz + WRITE(6,*)"TERMINATING PROGRAM!!!!" + STOP + ENDIF + IF(irx.eq.nnx)irx=irx-1 + IF(irz.eq.nnz)irz=irz-1 +! +! Location of receiver successfully found within the grid. Now approximate +! traveltime at receiver using bilinear interpolation from four +! surrounding grid points. Note that bilinear interpolation is a poor +! approximation when traveltime gradient varies significantly across a cell, +! particularly near the source. Thus, we use an improved approximation in this +! case. First, locate current source cell. +! + isx=INT((scx-gox)/dnx)+1 + isz=INT((scz-goz)/dnz)+1 + dpl=dnx*earth + rd1=dnz*earth*SIN(gox) + IF(rd1.LT.dpl)dpl=rd1 + rd1=dnz*earth*SIN(gox+(nnx-1)*dnx) + IF(rd1.LT.dpl)dpl=rd1 + sred=((scx-rcx1)*earth)**2 + sred=sred+((scz-rcz1)*earth*SIN(rcx1))**2 + sred=SQRT(sred) + IF(sred.LT.dpl)sw=1 + IF(isx.EQ.irx)THEN + IF(isz.EQ.irz)sw=1 + ENDIF + IF(sw.EQ.1)THEN +! +! Compute velocity at source and receiver +! + DO k=1,2 + DO l=1,2 + vss(k,l)=veln(isz-1+l,isx-1+k) + ENDDO + ENDDO + drx=(scx-gox)-(isx-1)*dnx + drz=(scz-goz)-(isz-1)*dnz + CALL bilinear(vss,drx,drz,vels) + DO k=1,2 + DO l=1,2 + vss(k,l)=veln(irz-1+l,irx-1+k) + ENDDO + ENDDO + drx=(rcx1-gox)-(irx-1)*dnx + drz=(rcz1-goz)-(irz-1)*dnz + CALL bilinear(vss,drx,drz,velr) + trr=2.0*sred/(vels+velr) + ELSE + drx=(rcx1-gox)-(irx-1)*dnx + drz=(rcz1-goz)-(irz-1)*dnz + trr=0.0 + DO k=1,2 + DO l=1,2 + produ=(1.0-ABS(((l-1)*dnz-drz)/dnz))*(1.0-ABS(((k-1)*dnx-drx)/dnx)) + trr=trr+ttn(irz-1+l,irx-1+k)*produ + ENDDO + ENDDO + ENDIF +! WRITE(10,*)yesray,trr +!!----------------------------------------------------------------- +! modified bu Hongjian Fang @ USTC +! count2=count2+1 +! cbst((no_p-1)*nsrc*nrc+(csid-1)*nrc+i)=trr + cbst1=trr +! call delsph(scx,scz,rcx(i),rcz(i),dist) +! travel_path(count2)=dist +!cbst_stat((no_p-1)*nsrc*nrc+(csid-1)*nrc+i)=yesray +!0605ENDDO +END SUBROUTINE srtimes + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: SUBROUTINE +! CODE: FORTRAN 90 +! This subroutine calculates ray path geometries for each +! source-receiver combination. It will also compute +! Frechet derivatives using these ray paths if required. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!SUBROUTINE rpaths(wrgf,csid,cfd,scx,scz) +!SUBROUTINE rpaths() +SUBROUTINE rpaths(scx,scz,fdm,surfrcx,surfrcz,writepath) +USE globalp +IMPLICIT NONE +INTEGER, PARAMETER :: i5=SELECTED_REAL_KIND(5,10) +INTEGER, PARAMETER :: nopath=0 +INTEGER :: i,j,k,l,m,n,ipx,ipz,ipxr,ipzr,nrp,sw +!fang!INTEGER :: wrgf,cfd,csid,ipxo,ipzo,isx,isz +INTEGER :: ipxo,ipzo,isx,isz +INTEGER :: ivx,ivz,ivxo,ivzo,nhp,maxrp +INTEGER :: ivxt,ivzt,ipxt,ipzt,isum,igref +INTEGER, DIMENSION (4) :: chp +REAL(KIND=i5) :: rayx,rayz +REAL(KIND=i10) :: dpl,rd1,rd2,xi,zi,vel,velo +REAL(KIND=i10) :: v,w,rigz,rigx,dinc,scx,scz +REAL(KIND=i10) :: dtx,dtz,drx,drz,produ,sred +REAL(KIND=i10), DIMENSION (:), ALLOCATABLE :: rgx,rgz +!fang!REAL(KIND=i5), DIMENSION (:,:), ALLOCATABLE :: fdm +REAL(KIND=i10), DIMENSION (4) :: vrat,vi,wi,vio,wio +!fang!------------------------------------------------ +real fdm(0:nvz+1,0:nvx+1) +REAL(KIND=i10) surfrcx,surfrcz +integer writepath +!fang!------------------------------------------------ +! +! ipx,ipz = Coordinates of cell containing current point +! ipxr,ipzr = Same as ipx,apz except for refined grid +! ipxo,ipzo = Coordinates of previous point +! rgx,rgz = (x,z) coordinates of ray geometry +! ivx,ivz = Coordinates of B-spline vertex containing current point +! ivxo,ivzo = Coordinates of previous point +! maxrp = maximum number of ray points +! nrp = number of points to describe ray +! dpl = incremental path length of ray +! xi,zi = edge of model coordinates +! dtx,dtz = components of gradT +! wrgf = Write out raypaths? (<0=all,0=no,>0=souce id) +! cfd = calculate Frechet derivatives? (0=no,1=yes) +! csid = current source id +! fdm = Frechet derivative matrix +! nhp = Number of ray segment-B-spline cell hit points +! vrat = length ratio of ray sub-segment +! chp = pointer to incremental change in x or z cell +! drx,drz = distance from reference node of cell +! produ = variable for trilinear interpolation +! vel = velocity at current point +! velo = velocity at previous point +! v,w = local variables of x,z +! vi,wi = B-spline basis functions at current point +! vio,wio = vi,wi for previous point +! ivxt,ivzt = temporary ivr,ivx,ivz values +! rigx,rigz = end point of sub-segment of ray path +! ipxt,ipzt = temporary ipx,ipz values +! dinc = path length of ray sub-segment +! rayr,rayx,rayz = ray path coordinates in single precision +! isx,isz = current source cell location +! scx,scz = current source coordinates +! sred = source to ray endpoint distance +! igref = ray endpoint lies in refined grid? (0=no,1=yes) +! nopath = switch to indicate that no path is present +! +! Allocate memory to arrays for storing ray path geometry +! +maxrp=nnx*nnz +ALLOCATE(rgx(maxrp+1), STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with ALLOCATE: SUBROUTINE rpaths: REAL rgx' +ENDIF +ALLOCATE(rgz(maxrp+1), STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with ALLOCATE: SUBROUTINE rpaths: REAL rgz' +ENDIF +! +! Allocate memory to partial derivative array +! +!fang!IF(cfd.EQ.1)THEN +!fang! ALLOCATE(fdm(0:nvz+1,0:nvx+1), STAT=checkstat) +!fang! IF(checkstat > 0)THEN +!fang! WRITE(6,*)'Error with ALLOCATE: SUBROUTINE rpaths: REAL fdm' +!fang! ENDIF +!fang!ENDIF +! +! Locate current source cell +! +IF(asgr.EQ.1)THEN + isx=INT((scx-goxr)/dnxr)+1 + isz=INT((scz-gozr)/dnzr)+1 +ELSE + isx=INT((scx-gox)/dnx)+1 + isz=INT((scz-goz)/dnz)+1 +ENDIF +! +! Set ray incremental path length equal to half width +! of cell +! + dpl=dnx*earth + rd1=dnz*earth*SIN(gox) + IF(rd1.LT.dpl)dpl=rd1 + rd1=dnz*earth*SIN(gox+(nnx-1)*dnx) + IF(rd1.LT.dpl)dpl=rd1 + dpl=0.5*dpl +! +! Loop through all the receivers +! +!fang!DO i=1,nrc +! +! If path does not exist, then cycle the loop +! +fdm=0 +!fang! IF(cfd.EQ.1)THEN +!fang! fdm=0.0 +!fang! ENDIF +!fang! IF(srs(i,csid).EQ.0)THEN +!fang! IF(wrgf.EQ.csid.OR.wrgf.LT.0)THEN +!fang! WRITE(40)nopath +!fang! ENDIF +!fang! IF(cfd.EQ.1)THEN +!fang! WRITE(50)nopath +!fang! ENDIF +!fang! CYCLE +!fang! ENDIF +! +! The first step is to locate the receiver in the grid. +! + ipx=INT((surfrcx-gox)/dnx)+1 + ipz=INT((surfrcz-goz)/dnz)+1 + sw=0 + IF(ipx.lt.1.or.ipx.ge.nnx)sw=1 + IF(ipz.lt.1.or.ipz.ge.nnz)sw=1 + IF(sw.eq.1)then + ipx=90.0-ipx*180.0/pi + ipz=ipz*180.0/pi + WRITE(6,*)"rpath Receiver lies outside model (lat,long)= ",ipx,ipz + WRITE(6,*)"TERMINATING PROGRAM!!!" + STOP + ENDIF + IF(ipx.eq.nnx)ipx=ipx-1 + IF(ipz.eq.nnz)ipz=ipz-1 +! +! First point of the ray path is the receiver +! + rgx(1)=surfrcx + rgz(1)=surfrcz +! +! Test to see if receiver is in source neighbourhood +! + sred=((scx-rgx(1))*earth)**2 + sred=sred+((scz-rgz(1))*earth*SIN(rgx(1)))**2 + sred=SQRT(sred) + IF(sred.LT.2.0*dpl)THEN + rgx(2)=scx + rgz(2)=scz + nrp=2 + sw=1 + ENDIF +! +! If required, see if receiver lies within refined grid +! + IF(asgr.EQ.1)THEN + ipxr=INT((surfrcx-goxr)/dnxr)+1 + ipzr=INT((surfrcz-gozr)/dnzr)+1 + igref=1 + IF(ipxr.LT.1.OR.ipxr.GE.nnxr)igref=0 + IF(ipzr.LT.1.OR.ipzr.GE.nnzr)igref=0 + IF(igref.EQ.1)THEN + IF(nstsr(ipzr,ipxr).NE.0.OR.nstsr(ipzr+1,ipxr).NE.0)igref=0 + IF(nstsr(ipzr,ipxr+1).NE.0.OR.nstsr(ipzr+1,ipxr+1).NE.0)igref=0 + ENDIF + ELSE + igref=0 + ENDIF +! +! Due to the method for calculating traveltime gradient, if the +! the ray end point lies in the source cell, then we are also done. +! + IF(sw.EQ.0)THEN + IF(asgr.EQ.1)THEN + IF(igref.EQ.1)THEN + IF(ipxr.EQ.isx)THEN + IF(ipzr.EQ.isz)THEN + rgx(2)=scx + rgz(2)=scz + nrp=2 + sw=1 + ENDIF + ENDIF + ENDIF + ELSE + IF(ipx.EQ.isx)THEN + IF(ipz.EQ.isz)THEN + rgx(2)=scx + rgz(2)=scz + nrp=2 + sw=1 + ENDIF + ENDIF + ENDIF + ENDIF +! +! Now trace ray from receiver to "source" +! + DO j=1,maxrp + IF(sw.EQ.1)EXIT +! +! Calculate traveltime gradient vector for current cell using +! a first-order or second-order scheme. +! + IF(igref.EQ.1)THEN +! +! In this case, we are in the refined grid. +! +! First order scheme applied here. +! + dtx=ttnr(ipzr,ipxr+1)-ttnr(ipzr,ipxr) + dtx=dtx+ttnr(ipzr+1,ipxr+1)-ttnr(ipzr+1,ipxr) + dtx=dtx/(2.0*earth*dnxr) + dtz=ttnr(ipzr+1,ipxr)-ttnr(ipzr,ipxr) + dtz=dtz+ttnr(ipzr+1,ipxr+1)-ttnr(ipzr,ipxr+1) + dtz=dtz/(2.0*earth*SIN(rgx(j))*dnzr) + ELSE +! +! Here, we are in the coarse grid. +! +! First order scheme applied here. +! + dtx=ttn(ipz,ipx+1)-ttn(ipz,ipx) + dtx=dtx+ttn(ipz+1,ipx+1)-ttn(ipz+1,ipx) + dtx=dtx/(2.0*earth*dnx) + dtz=ttn(ipz+1,ipx)-ttn(ipz,ipx) + dtz=dtz+ttn(ipz+1,ipx+1)-ttn(ipz,ipx+1) + dtz=dtz/(2.0*earth*SIN(rgx(j))*dnz) + ENDIF +! +! Calculate the next ray path point +! + rd1=SQRT(dtx**2+dtz**2) + rgx(j+1)=rgx(j)-dpl*dtx/(earth*rd1) + rgz(j+1)=rgz(j)-dpl*dtz/(earth*SIN(rgx(j))*rd1) +! +! Determine which cell the new ray endpoint +! lies in. +! + ipxo=ipx + ipzo=ipz + IF(asgr.EQ.1)THEN +! +! Here, we test to see whether the ray endpoint lies +! within a cell of the refined grid +! + ipxr=INT((rgx(j+1)-goxr)/dnxr)+1 + ipzr=INT((rgz(j+1)-gozr)/dnzr)+1 + igref=1 + IF(ipxr.LT.1.OR.ipxr.GE.nnxr)igref=0 + IF(ipzr.LT.1.OR.ipzr.GE.nnzr)igref=0 + IF(igref.EQ.1)THEN + IF(nstsr(ipzr,ipxr).NE.0.OR.nstsr(ipzr+1,ipxr).NE.0)igref=0 + IF(nstsr(ipzr,ipxr+1).NE.0.OR.nstsr(ipzr+1,ipxr+1).NE.0)igref=0 + ENDIF + ipx=INT((rgx(j+1)-gox)/dnx)+1 + ipz=INT((rgz(j+1)-goz)/dnz)+1 + ELSE + ipx=INT((rgx(j+1)-gox)/dnx)+1 + ipz=INT((rgz(j+1)-goz)/dnz)+1 + igref=0 + ENDIF +! +! Test the proximity of the source to the ray end point. +! If it is less than dpl then we are done +! + sred=((scx-rgx(j+1))*earth)**2 + sred=sred+((scz-rgz(j+1))*earth*SIN(rgx(j+1)))**2 + sred=SQRT(sred) + sw=0 + IF(sred.LT.2.0*dpl)THEN + rgx(j+2)=scx + rgz(j+2)=scz + nrp=j+2 + sw=1 +!fang! IF(cfd.NE.1)EXIT + ENDIF +! +! Due to the method for calculating traveltime gradient, if the +! the ray end point lies in the source cell, then we are also done. +! + IF(sw.EQ.0)THEN + IF(asgr.EQ.1)THEN + IF(igref.EQ.1)THEN + IF(ipxr.EQ.isx)THEN + IF(ipzr.EQ.isz)THEN + rgx(j+2)=scx + rgz(j+2)=scz + nrp=j+2 + sw=1 + !fang! IF(cfd.NE.1)EXIT + ENDIF + ENDIF + ENDIF + ELSE + IF(ipx.EQ.isx)THEN + IF(ipz.EQ.isz)THEN + rgx(j+2)=scx + rgz(j+2)=scz + nrp=j+2 + sw=1 + !fang! IF(cfd.NE.1)EXIT + ENDIF + ENDIF + ENDIF + ENDIF +! +! Test whether ray path segment extends beyond +! box boundaries +! + IF(ipx.LT.1)THEN + rgx(j+1)=gox + ipx=1 + rbint=1 + ENDIF + IF(ipx.GE.nnx)THEN + rgx(j+1)=gox+(nnx-1)*dnx + ipx=nnx-1 + rbint=1 + ENDIF + IF(ipz.LT.1)THEN + rgz(j+1)=goz + ipz=1 + rbint=1 + ENDIF + IF(ipz.GE.nnz)THEN + rgz(j+1)=goz+(nnz-1)*dnz + ipz=nnz-1 + rbint=1 + ENDIF +! +! Calculate the Frechet derivatives if required. +! + !fang! IF(cfd.EQ.1)THEN +! +! First determine which B-spline cell the refined cells +! containing the ray path segment lies in. If they lie +! in more than one, then we need to divide the problem +! into separate parts (up to three). +! + ivx=INT((ipx-1)/gdx)+1 + ivz=INT((ipz-1)/gdz)+1 + ivxo=INT((ipxo-1)/gdx)+1 + ivzo=INT((ipzo-1)/gdz)+1 +! +! Calculate up to two hit points between straight +! ray segment and cell faces. +! + nhp=0 + IF(ivx.NE.ivxo)THEN + nhp=nhp+1 + IF(ivx.GT.ivxo)THEN + xi=gox+(ivx-1)*dvx + ELSE + xi=gox+ivx*dvx + ENDIF + vrat(nhp)=(xi-rgx(j))/(rgx(j+1)-rgx(j)) + chp(nhp)=1 + ENDIF + IF(ivz.NE.ivzo)THEN + nhp=nhp+1 + IF(ivz.GT.ivzo)THEN + zi=goz+(ivz-1)*dvz + ELSE + zi=goz+ivz*dvz + ENDIF + rd1=(zi-rgz(j))/(rgz(j+1)-rgz(j)) + IF(nhp.EQ.1)THEN + vrat(nhp)=rd1 + chp(nhp)=2 + ELSE + IF(rd1.GE.vrat(nhp-1))THEN + vrat(nhp)=rd1 + chp(nhp)=2 + ELSE + vrat(nhp)=vrat(nhp-1) + chp(nhp)=chp(nhp-1) + vrat(nhp-1)=rd1 + chp(nhp-1)=2 + ENDIF + ENDIF + ENDIF + nhp=nhp+1 + vrat(nhp)=1.0 + chp(nhp)=0 +! +! Calculate the velocity, v and w values of the +! first point +! + drx=(rgx(j)-gox)-(ipxo-1)*dnx + drz=(rgz(j)-goz)-(ipzo-1)*dnz + vel=0.0 + DO l=1,2 + DO m=1,2 + produ=(1.0-ABS(((m-1)*dnz-drz)/dnz)) + produ=produ*(1.0-ABS(((l-1)*dnx-drx)/dnx)) + IF(ipzo-1+m.LE.nnz.AND.ipxo-1+l.LE.nnx)THEN + vel=vel+veln(ipzo-1+m,ipxo-1+l)*produ + ENDIF + ENDDO + ENDDO + drx=(rgx(j)-gox)-(ivxo-1)*dvx + drz=(rgz(j)-goz)-(ivzo-1)*dvz + v=drx/dvx + w=drz/dvz +! +! Calculate the 12 basis values at the point +! + vi(1)=(1.0-v)**3/6.0 + vi(2)=(4.0-6.0*v**2+3.0*v**3)/6.0 + vi(3)=(1.0+3.0*v+3.0*v**2-3.0*v**3)/6.0 + vi(4)=v**3/6.0 + wi(1)=(1.0-w)**3/6.0 + wi(2)=(4.0-6.0*w**2+3.0*w**3)/6.0 + wi(3)=(1.0+3.0*w+3.0*w**2-3.0*w**3)/6.0 + wi(4)=w**3/6.0 + ivxt=ivxo + ivzt=ivzo +! +! Now loop through the one or more sub-segments of the +! ray path segment and calculate partial derivatives +! + DO k=1,nhp + velo=vel + vio=vi + wio=wi + IF(k.GT.1)THEN + IF(chp(k-1).EQ.1)THEN + ivxt=ivx + ELSE IF(chp(k-1).EQ.2)THEN + ivzt=ivz + ENDIF + ENDIF +! +! Calculate the velocity, v and w values of the +! new point +! + rigz=rgz(j)+vrat(k)*(rgz(j+1)-rgz(j)) + rigx=rgx(j)+vrat(k)*(rgx(j+1)-rgx(j)) + ipxt=INT((rigx-gox)/dnx)+1 + ipzt=INT((rigz-goz)/dnz)+1 + drx=(rigx-gox)-(ipxt-1)*dnx + drz=(rigz-goz)-(ipzt-1)*dnz + vel=0.0 + DO m=1,2 + DO n=1,2 + produ=(1.0-ABS(((n-1)*dnz-drz)/dnz)) + produ=produ*(1.0-ABS(((m-1)*dnx-drx)/dnx)) + IF(ipzt-1+n.LE.nnz.AND.ipxt-1+m.LE.nnx)THEN + vel=vel+veln(ipzt-1+n,ipxt-1+m)*produ + ENDIF + ENDDO + ENDDO + drx=(rigx-gox)-(ivxt-1)*dvx + drz=(rigz-goz)-(ivzt-1)*dvz + v=drx/dvx + w=drz/dvz +! +! Calculate the 8 basis values at the new point +! + vi(1)=(1.0-v)**3/6.0 + vi(2)=(4.0-6.0*v**2+3.0*v**3)/6.0 + vi(3)=(1.0+3.0*v+3.0*v**2-3.0*v**3)/6.0 + vi(4)=v**3/6.0 + wi(1)=(1.0-w)**3/6.0 + wi(2)=(4.0-6.0*w**2+3.0*w**3)/6.0 + wi(3)=(1.0+3.0*w+3.0*w**2-3.0*w**3)/6.0 + wi(4)=w**3/6.0 +! +! Calculate the incremental path length +! + IF(k.EQ.1)THEN + dinc=vrat(k)*dpl + ELSE + dinc=(vrat(k)-vrat(k-1))*dpl + ENDIF +! +! Now compute the 16 contributions to the partial +! derivatives. +! + DO l=1,4 + DO m=1,4 + rd1=vi(m)*wi(l)/vel**2 + rd2=vio(m)*wio(l)/velo**2 + rd1=-(rd1+rd2)*dinc/2.0 + !fang! rd1=vi(m)*wi(l) + !fang! rd2=vio(m)*wio(l) + !fang! rd1=(rd1+rd2)*dinc/2.0 + rd2=fdm(ivzt-2+l,ivxt-2+m) + fdm(ivzt-2+l,ivxt-2+m)=rd1+rd2 + ENDDO + ENDDO + ENDDO + !fang! ENDIF +!fang! IF(j.EQ.maxrp.AND.sw.EQ.0)THEN +!fang! WRITE(6,*)'Error with ray path detected!!!' +!fang! WRITE(6,*)'Source id: ',csid +!fang! WRITE(6,*)'Receiver id: ',i +!fang! ENDIF + ENDDO +! +! Write ray paths to output file +! +!fang! IF(wrgf.EQ.csid.OR.wrgf.LT.0)THEN + if(writepath == 1) then + WRITE(40,*)'#',nrp + DO j=1,nrp + rayx=(pi/2-rgx(j))*180.0/pi + rayz=rgz(j)*180.0/pi + WRITE(40,*)rayx,rayz + ENDDO + endif +!fang! ENDIF +! +! Write partial derivatives to output file +! +!fang! IF(cfd.EQ.1)THEN +!fang!! +!fang!! Determine the number of non-zero elements. +!fang!! +!fang! isum=0 +!fang! DO j=0,nvz+1 +!fang! DO k=0,nvx+1 +!fang! IF(ABS(fdm(j,k)).GE.ftol)isum=isum+1 +!fang! ENDDO +!fang! ENDDO +!fang! WRITE(50)isum +!fang! isum=0 +!fang! DO j=0,nvz+1 +!fang! DO k=0,nvx+1 +!fang! isum=isum+1 +!fang! IF(ABS(fdm(j,k)).GE.ftol)WRITE(50)isum,fdm(j,k) +!fang! ENDDO +!fang! ENDDO +!fang! ENDIF +!fang!ENDDO +!fang!IF(cfd.EQ.1)THEN +!fang! DEALLOCATE(fdm, STAT=checkstat) +!fang! IF(checkstat > 0)THEN +!fang! WRITE(6,*)'Error with DEALLOCATE: SUBROUTINE rpaths: fdm' +!fang! ENDIF +!fang!ENDIF +DEALLOCATE(rgx,rgz, STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with DEALLOCATE: SUBROUTINE rpaths: rgx,rgz' +ENDIF +END SUBROUTINE rpaths + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! TYPE: SUBROUTINE +! CODE: FORTRAN 90 +! This subroutine is passed four node values which lie on +! the corners of a rectangle and the coordinates of a point +! lying within the rectangle. It calculates the value at +! the internal point by using bilinear interpolation. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +SUBROUTINE bilinear(nv,dsx,dsz,biv) +USE globalp +IMPLICIT NONE +INTEGER :: i,j +REAL(KIND=i10) :: dsx,dsz,biv +REAL(KIND=i10), DIMENSION(2,2) :: nv +REAL(KIND=i10) :: produ +! +! nv = four node vertex values +! dsx,dsz = distance between internal point and top left node +! dnx,dnz = width and height of node rectangle +! biv = value at internal point calculated by bilinear interpolation +! produ = product variable +! +biv=0.0 +DO i=1,2 + DO j=1,2 + produ=(1.0-ABS(((i-1)*dnx-dsx)/dnx))*(1.0-ABS(((j-1)*dnz-dsz)/dnz)) + biv=biv+nv(i,j)*produ + ENDDO +ENDDO +END SUBROUTINE bilinear + + + subroutine refineGrid2LayerMdl(minthk0,mmax,dep,vp,vs,rho,& + rmax,rdep,rvp,rvs,rrho,rthk) +!!--------------------------------------------------------------------c +!!refine grid based model to layerd based model +!!INPUT: minthk: is the minimum thickness of the refined layered model +!! mmax: number of depth grid points in the model +!! dep, vp, vs, rho: the depth-grid model parameters +!! rmax: number of layers in the fined layered model +!! rdep, rvp, rvs, rrho, rthk: the refined layered velocity model +!! + implicit none + integer NL + parameter (NL=200) + integer mmax,rmax + real minthk0 + real minthk + real dep(*),vp(*),vs(*),rho(*) + real rdep(NL),rvp(NL),rvs(NL),rrho(NL),rthk(NL) + integer nsublay(NL) + real thk,newthk,initdep + integer i,j,k,ngrid + + k = 0 + initdep = 0.0 + do i = 1, mmax-1 + thk = dep(i+1)-dep(i) + minthk = thk/minthk0 + nsublay(i) = int((thk+1.0e-4)/minthk) + 1 + ngrid = nsublay(i)+1 + newthk = thk/nsublay(i) + do j = 1, nsublay(i) + k = k + 1 + rthk(k) = newthk + rdep(k) = initdep + rthk(k) + initdep = rdep(k) + rvp(k) = vp(i)+(2*j-1)*(vp(i+1)-vp(i))/(2*nsublay(i)) + rvs(k) = vs(i)+(2*j-1)*(vs(i+1)-vs(i))/(2*nsublay(i)) + rrho(k) = rho(i)+(2*j-1)*(rho(i+1)-rho(i))/(2*nsublay(i)) + enddo + enddo +!! half space model + k = k + 1 + rthk(k) = 0.0 + rvp(k) = vp(mmax) + rvs(k) = vs(mmax) + rrho(k) = rho(mmax) + + rmax = k + +!! do i = 1, mmax +!! write(*,*) dep(i),vp(i),vs(i),rho(i) +!! enddo +!! print *, '---------------------------------' +!! do i = 1, rmax +!! write(*,*) rdep(i),rthk(i),rvp(i),rvs(i),rrho(i) +!! enddo + + return + end +subroutine synthetic(nx,ny,nz,nparpi,vels,obst, & + goxdf,gozdf,dvxdf,dvzdf,kmaxRc,kmaxRg,kmaxLc,kmaxLg, & + tRc,tRg,tLc,tLg,wavetype,igrt,periods,depz,minthk, & + scxf,sczf,rcxf,rczf,nrc1,nsrcsurf1,knum1,kmax,nsrcsurf,nrcf,noiselevel) + +USE globalp +USE traveltime +IMPLICIT NONE +!CHARACTER (LEN=30) ::grid,frechet +!CHARACTER (LEN=40) :: sources,receivers,otimes +!CHARACTER (LEN=30) :: travelt,rtravel,wrays,cdum +INTEGER :: i,j,k,l,nsrc,tnr,urg +INTEGER :: sgs,isx,isz,sw,idm1,idm2,nnxb,nnzb +INTEGER :: ogx,ogz,grdfx,grdfz,maxbt +REAL(KIND=i10) :: x,z,goxb,gozb,dnxb,dnzb +!REAL(KIND=i10), DIMENSION (:,:), ALLOCATABLE :: scxf,sczf +!REAL(KIND=i10), DIMENSION (:,:,:), ALLOCATABLE :: rcxf,rczf +! +! sources = File containing source locations +! receivers = File containing receiver locations +! grid = File containing grid of velocity vertices for +! resampling on a finer grid with cubic B-splines +! frechet = output file containing matrix of frechet derivatives +! travelt = File name for storage of traveltime field +! wttf = Write traveltimes to file? (0=no,>0=source id) +! fom = Use first-order(0) or mixed-order(1) scheme +! nsrc = number of sources +! scx,scz = source location in r,x,z +! scx,scz = source location in r,x,z +! x,z = temporary variables for source location +! fsrt = find source-receiver traveltimes? (0=no,1=yes) +! rtravel = output file for source-receiver traveltimes +! cdum = dummy character variable ! wrgf = write ray geometries to file? (<0=all,0=no,>0=source id.) +! wrays = file containing raypath geometries +! cfd = calculate Frechet derivatives? (0=no, 1=yes) +! tnr = total number of receivers +! sgs = Extent of refined source grid +! isx,isz = cell containing source +! nnxb,nnzb = Backup for nnz,nnx +! goxb,gozb = Backup for gox,goz +! dnxb,dnzb = Backup for dnx,dnz +! ogx,ogz = Location of refined grid origin +! gridfx,grdfz = Number of refined nodes per cell +! urg = use refined grid (0=no,1=yes,2=previously used) +! maxbt = maximum size of narrow band binary tree +! otimes = file containing source-receiver association information +!c----------------------------------------------------------------- +! variables defined by Hongjian Fang + integer nx,ny,nz + integer kmax,nsrcsurf,nrcf + real vels(nx,ny,nz) + real obst(*) + real goxdf,gozdf,dvxdf,dvzdf + integer kmaxRc,kmaxRg,kmaxLc,kmaxLg + real*8 tRc(*),tRg(*),tLc(*),tLg(*) + integer wavetype(nsrcsurf,kmax) + integer periods(nsrcsurf,kmax),nrc1(nsrcsurf,kmax),nsrcsurf1(kmax) + integer knum1(kmax),igrt(nsrcsurf,kmax) + real scxf(nsrcsurf,kmax),sczf(nsrcsurf,kmax),rcxf(nrcf,nsrcsurf,kmax),rczf(nrcf,nsrcsurf,kmax) + real minthk + integer nparpi + + + real vpz(nz),vsz(nz),rhoz(nz),depz(nz) + real*8 pvRc(nx*ny,kmaxRc),pvRg(nx*ny,kmaxRg),pvLc(nx*ny,kmaxLc),pvLg(nx*ny,kmaxLg) + real*8 velf(ny*nx) + integer kmax1,kmax2,kmax3,count1 + integer igr + integer iwave + integer knumi,srcnum + real cbst1 + real noiselevel + real gaussian + external gaussian + integer ii,jj,kk,nn,istep +gdx=5 +gdz=5 +asgr=1 +sgdl=8 +sgs=8 +earth=6371.0 +fom=1 +snb=0.5 +goxd=goxdf +gozd=gozdf +dvxd=dvxdf +dvzd=dvzdf +nvx=nx-2 +nvz=ny-2 +ALLOCATE(velv(0:nvz+1,0:nvx+1), STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with ALLOCATE: SUBROUTINE gridder: REAL velv' +ENDIF +! +! Convert from degrees to radians +! +dvx=dvxd*pi/180.0 +dvz=dvzd*pi/180.0 +gox=(90.0-goxd)*pi/180.0 +goz=gozd*pi/180.0 +! +! Compute corresponding values for propagation grid. +! +nnx=(nvx-1)*gdx+1 +nnz=(nvz-1)*gdz+1 +dnx=dvx/gdx +dnz=dvz/gdz +dnxd=dvxd/gdx +dnzd=dvzd/gdz +ALLOCATE(veln(nnz,nnx), STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with ALLOCATE: SUBROUTINE gridder: REAL veln' +ENDIF + +! +! Call a subroutine which reads in the velocity grid +! +!CALL gridder(grid) +! +! Read in all source coordinates. +! +! +! Now work out, source by source, the first-arrival traveltime +! field plus source-receiver traveltimes +! and ray paths if required. First, allocate memory to the +! traveltime field array +! +ALLOCATE(ttn(nnz,nnx), STAT=checkstat) +IF(checkstat > 0)THEN + WRITE(6,*)'Error with ALLOCATE: PROGRAM fmmin2d: REAL ttn' +ENDIF + rbint=0 +! +! Allocate memory for node status and binary trees +! +ALLOCATE(nsts(nnz,nnx)) +maxbt=NINT(snb*nnx*nnz) +ALLOCATE(btg(maxbt)) + +!allocate(fdm(0:nvz+1,0:nvx+1)) + + if(kmaxRc.gt.0) then + iwave=2 + igr=0 + call caldespersion(nx,ny,nz,vels,pvRc, & + iwave,igr,kmaxRc,tRc,depz,minthk) + endif + + if(kmaxRg.gt.0) then + iwave=2 + igr=1 + call caldespersion(nx,ny,nz,vels,pvRg, & + iwave,igr,kmaxRg,tRg,depz,minthk) + endif + + if(kmaxLc.gt.0) then + iwave=1 + igr=0 + call caldespersion(nx,ny,nz,vels,pvLc, & + iwave,igr,kmaxLc,tLc,depz,minthk) + endif + + if(kmaxLg.gt.0) then + iwave=1 + igr=1 + call caldespersion(nx,ny,nz,vels,pvLg, & + iwave,igr,kmaxLg,tLg,depz,minthk) + endif + +!nar=0 +count1=0 + +!sen_vs=0 +!sen_vp=0 +!sen_rho=0 +kmax1=kmaxRc +kmax2=kmaxRc+kmaxRg +kmax3=kmaxRc+kmaxRg+kmaxLc +do knumi=1,kmax +do srcnum=1,nsrcsurf1(knum1(knumi)) + if(wavetype(srcnum,knum1(knumi))==2.and.igrt(srcnum,knum1(knumi))==0) then + velf(1:nx*ny)=pvRc(1:nx*ny,periods(srcnum,knum1(knumi))) +! sen_vs(:,1:kmax1,:)=sen_vsRc(:,1:kmaxRc,:)!(:,nt(istep),:) +! sen_vp(:,1:kmax1,:)=sen_vpRc(:,1:kmaxRc,:)!(:,nt(istep),:) +! sen_rho(:,1:kmax1,:)=sen_rhoRc(:,1:kmaxRc,:)!(:,nt(istep),:) + endif + if(wavetype(srcnum,knum1(knumi))==2.and.igrt(srcnum,knum1(knumi))==1) then + velf(1:nx*ny)=pvRg(1:nx*ny,periods(srcnum,knum1(knumi))) +! sen_vs(:,kmax1+1:kmax2,:)=sen_vsRg(:,1:kmaxRg,:)!(:,nt,:) +! sen_vp(:,kmax1+1:kmax2,:)=sen_vpRg(:,1:kmaxRg,:)!(:,nt,:) +! sen_rho(:,kmax1+1:kmax2,:)=sen_rhoRg(:,1:kmaxRg,:)!(:,nt,:) + endif + if(wavetype(srcnum,knum1(knumi))==1.and.igrt(srcnum,knum1(knumi))==0) then + velf(1:nx*ny)=pvLc(1:nx*ny,periods(srcnum,knum1(knumi))) +! sen_vs(:,kmax2+1:kmax3,:)=sen_vsLc(:,1:kmaxLc,:)!(:,nt,:) +! sen_vp(:,kmax2+1:kmax3,:)=sen_vpLc(:,1:kmaxLc,:)!(:,nt,:) +! sen_rho(:,kmax2+1:kmax3,:)=sen_rhoLc(:,1:kmaxLc,:)!(:,nt,:) + endif + if(wavetype(srcnum,knum1(knumi))==1.and.igrt(srcnum,knum1(knumi))==1) then + velf(1:nx*ny)=pvLg(1:nx*ny,periods(srcnum,knum1(knumi))) +! sen_vs(:,kmax3+1:kmax,:)=sen_vsLg(:,1:kmaxLg,:)!(:,nt,:) +! sen_vp(:,kmax3+1:kmax,:)=sen_vpLg(:,1:kmaxLg,:)!(:,nt,:) +! sen_rho(:,kmax3+1:kmax,:)=sen_rhoLg(:,1:kmaxLg,:)!(:,nt,:) + endif + +call gridder(velf) + x=scxf(srcnum,knum1(knumi)) + z=sczf(srcnum,knum1(knumi)) +! +! Begin by computing refined source grid if required +! + urg=0 + IF(asgr.EQ.1)THEN +! +! Back up coarse velocity grid to a holding matrix +! + ALLOCATE(velnb(nnz,nnx)) +! MODIFIEDY BY HONGJIAN FANG @ USTC 2014/04/17 + velnb(1:nnz,1:nnx)=veln(1:nnz,1:nnx) + nnxb=nnx + nnzb=nnz + dnxb=dnx + dnzb=dnz + goxb=gox + gozb=goz +! +! Identify nearest neighbouring node to source +! + isx=INT((x-gox)/dnx)+1 + isz=INT((z-goz)/dnz)+1 + sw=0 + IF(isx.lt.1.or.isx.gt.nnx)sw=1 + IF(isz.lt.1.or.isz.gt.nnz)sw=1 + IF(sw.eq.1)then + isx=90.0-isx*180.0/pi + isz=isz*180.0/pi + WRITE(6,*)"Source lies outside bounds of model (lat,long)= ",isx,isz + WRITE(6,*)"TERMINATING PROGRAM!!!" + STOP + ENDIF + IF(isx.eq.nnx)isx=isx-1 + IF(isz.eq.nnz)isz=isz-1 +! +! Now find rectangular box that extends outward from the nearest source node +! to "sgs" nodes away. +! + vnl=isx-sgs + IF(vnl.lt.1)vnl=1 + vnr=isx+sgs + IF(vnr.gt.nnx)vnr=nnx + vnt=isz-sgs + IF(vnt.lt.1)vnt=1 + vnb=isz+sgs + IF(vnb.gt.nnz)vnb=nnz + nrnx=(vnr-vnl)*sgdl+1 + nrnz=(vnb-vnt)*sgdl+1 + drnx=dvx/REAL(gdx*sgdl) + drnz=dvz/REAL(gdz*sgdl) + gorx=gox+dnx*(vnl-1) + gorz=goz+dnz*(vnt-1) + nnx=nrnx + nnz=nrnz + dnx=drnx + dnz=drnz + gox=gorx + goz=gorz +! +! Reallocate velocity and traveltime arrays if nnx>nnxb or +! nnz 0)THEN + WRITE(6,*)'Error with DEALLOCATE: PROGRAM fmmin2d: velnb' + ENDIF +ENDIF +enddo +enddo +!deallocate(fdm) +deallocate(velv,veln,ttn,nsts,btg) +END subroutine +subroutine caldespersion(nx,ny,nz,vel,pvRc, & + iwave,igr,kmaxRc,tRc,depz,minthk) + use omp_lib + implicit none + + integer nx,ny,nz + real vel(nx,ny,nz) + + integer iwave,igr + real minthk + real depz(nz) + integer kmaxRc + real*8 tRc(kmaxRc) + real*8 pvRc(nx*ny,kmaxRc) + + + + real vpz(nz),vsz(nz),rhoz(nz) + integer mmax,iflsph,mode,rmax + integer ii,jj,k,i,nn,kk + integer,parameter::NL=200 + integer,parameter::NP=60 + real*8 cg1(NP),cg2(NP),cga,cgRc(NP) + real rdep(NL),rvp(NL),rvs(NL),rrho(NL),rthk(NL) + real depm(NL),vpm(NL),vsm(NL),rhom(NL),thkm(NL) + real dlnVs,dlnVp,dlnrho + + + mmax=nz + iflsph=1 + mode=1 + dlnVs=0.01 + dlnVp=0.01 + dlnrho=0.01 + +!$omp parallel & +!$omp default(private) & +!$omp shared(depz,nx,ny,nz,minthk,kmaxRc,mmax,vel) & +!$omp shared(tRc,pvRc,iflsph,iwave,mode,igr) +!$omp do + do jj=1,ny + do ii=1,nx + vsz(1:nz)=vel(ii,jj,1:nz) + ! some other emperical relationship maybe better, + do k=1,nz + vpz(k)=0.9409 + 2.0947*vsz(k) - 0.8206*vsz(k)**2+ & + 0.2683*vsz(k)**3 - 0.0251*vsz(k)**4 + rhoz(k)=1.6612*vpz(k) - 0.4721*vpz(k)**2 + & + 0.0671*vpz(k)**3 - 0.0043*vpz(k)**4 + & + 0.000106*vpz(k)**5 + enddo + + call refineGrid2LayerMdl(minthk,mmax,depz,vpz,vsz,rhoz,rmax,rdep,& + rvp,rvs,rrho,rthk) + call surfdisp96(rthk,rvp,rvs,rrho,rmax,iflsph,iwave,mode,igr,kmaxRc,& + tRc,cgRc) + pvRc((jj-1)*nx+ii,1:kmaxRc)=cgRc(1:kmaxRc) + enddo + enddo +!$omp end do +!$omp end parallel + end subroutine + + diff --git a/src/DSurfTomo b/src/DSurfTomo new file mode 100755 index 0000000..74536d5 Binary files /dev/null and b/src/DSurfTomo differ diff --git a/src/Makefile b/src/Makefile new file mode 100644 index 0000000..06cdfb2 --- /dev/null +++ b/src/Makefile @@ -0,0 +1,20 @@ +CMD = DSurfTomo +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\ + aprod.f90 gaussian.f90 main.f90 +FSRCS = surfdisp96.f +OBJS = $(F90SRCS:%.f90=%.o) $(FSRCS:%.f=%.o) CalSurfG.o +all:$(CMD) +$(CMD):$(OBJS) + $(FC) -fopenmp $^ -o $@ +CalSurfG.o:CalSurfG.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) diff --git a/src/aprod.f90 b/src/aprod.f90 new file mode 100644 index 0000000..5e17045 --- /dev/null +++ b/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 diff --git a/src/delsph.f90 b/src/delsph.f90 new file mode 100644 index 0000000..c9f7170 --- /dev/null +++ b/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 diff --git a/src/gaussian.f90 b/src/gaussian.f90 new file mode 100644 index 0000000..4cb5775 --- /dev/null +++ b/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 diff --git a/src/lsmrDataModule.f90 b/src/lsmrDataModule.f90 new file mode 100644 index 0000000..fb94f29 --- /dev/null +++ b/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 diff --git a/src/lsmrModule.f90 b/src/lsmrModule.f90 new file mode 100644 index 0000000..395ef00 --- /dev/null +++ b/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 +! Michael Saunders +! 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 diff --git a/src/lsmrblas.f90 b/src/lsmrblas.f90 new file mode 100644 index 0000000..31574e2 --- /dev/null +++ b/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 +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/src/lsmrblasInterface.f90 b/src/lsmrblasInterface.f90 new file mode 100644 index 0000000..58cefa0 --- /dev/null +++ b/src/lsmrblasInterface.f90 @@ -0,0 +1,41 @@ +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! File lsmrblasInterface.f90 +! +! BLAS1 Interfaces: ddot dnrm2 dscal +! +! Maintained by Michael Saunders . +! +! 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 diff --git a/src/main.f90 b/src/main.f90 new file mode 100644 index 0000000..1500f88 --- /dev/null +++ b/src/main.f90 @@ -0,0 +1,616 @@ + ! 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 + 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 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 + 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,*) maxiter + 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,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 + + 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 + + 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' + + + 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/sqrt(real(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/sqrt(real(dall)) + + 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 diff --git a/src/surfdisp96.f b/src/surfdisp96.f new file mode 100644 index 0000000..1e61103 --- /dev/null +++ b/src/surfdisp96.f @@ -0,0 +1,1062 @@ +c----------------------------------------------------------------------c +c c +c COMPUTER PROGRAMS IN SEISMOLOGY c +c VOLUME IV c +c c +c PROGRAM: SRFDIS c +c c +c COPYRIGHT 1986, 1991 c +c D. R. Russell, R. B. Herrmann c +c Department of Earth and Atmospheric Sciences c +c Saint Louis University c +c 221 North Grand Boulevard c +c St. Louis, Missouri 63103 c +c U. S. A. c +c c +c----------------------------------------------------------------------c +c This is a combination of program 'surface80' which search the poles +c on C-T domain, and the program 'surface81' which search in the F-K +c domain. The input data is slightly different with its precessors. +c -Wang 06/06/83. +c +c The program calculates the dispersion values for any +c layered model, any frequency, and any mode. +c +c This program will accept one liquid layer at the surface. +c In such case ellipticity of rayleigh wave is that at the +c top of solid array. Love wave communications ignore +c liquid layer. +c +c Program developed by Robert B Herrmann Saint Louis +c univ. Nov 1971, and revised by C. Y. Wang on Oct 1981. +c Modified for use in surface wave inversion, and +c addition of spherical earth flattening transformation, by +c David R. Russell, St. Louis University, Jan. 1984. +c +c Changes +c 28 JAN 2003 - fixed minor but for sphericity correction by +c saving one parameter in subroutine sphere +c 20 JUL 2004 - removed extraneous line at line 550 +c since dc not defined +c if(dabs(c1-c2) .le. dmin1(1.d-6*c1,0.005d+0*dc) )go to 1000 +c 28 DEC 2007 - changed the Earth flattening to now use layer +c midpoint and the Biswas (1972: PAGEOPH 96, 61-74, 1972) +c density mapping for P-SV - note a true comparison +c requires the ability to handle a fluid core for SH and SV +c Also permit one layer with fluid is base of the velocity is 0.001 km/sec +c----- +c 13 JAN 2010 - modified by Huajian Yao at MIT for calculation of +c group or phase velocities +c----- + + subroutine surfdisp96(thkm,vpm,vsm,rhom,nlayer,iflsph,iwave, + & mode,igr,kmax,t,cg) + + parameter(LER=0,LIN=5,LOT=66) + integer NL, NL2, NLAY + parameter(NL=200,NLAY=200,NL2=NL+NL) + integer NP + parameter (NP=60) + +c----- +c LIN - unit for FORTRAN read from terminal +c LOT - unit for FORTRAN write to terminal +c LER - unit for FORTRAN error output to terminal +c NL - layers in model +c NP - number of unique periods +c----- +c----- parameters +c thkm, vpm, vsm, rhom: model for dispersion calculation +c nlayer - I4: number of layers in the model +c iflsph - I4: 0 flat earth model, 1 spherical earth model +c iwave - I4: 1 Love wave, 2 Rayleigh wave +c mode - I4: ith mode of surface wave, 1 fundamental, 2 first higher, .... +c igr - I4: 0 phase velocity, > 0 group velocity +c kmax - I4: number of periods (t) for dispersion calculation +c t - period vector (t(NP)) +c cg - output phase or group velocities (vector,cg(NP)) +c----- + real*4 thkm(NLAY),vpm(NLAY),vsm(NLAY),rhom(NLAY) + integer nlayer,iflsph,iwave,mode,igr,kmax + double precision twopi,one,onea + double precision cc,c1,clow,cm,dc,t1 + double precision t(NP),c(NP),cb(NP),cg(NP) + real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL) +c common/modl/ d,a,b,rho,rtp,dtp,btp +c common/para/ mmax,llw,twopi + integer*4 iverb(2) + integer*4 llw + integer*4 nsph, ifunc, idispl, idispr, is, ie + real*4 sone0, ddc0, h0, sone, ddc, h + +c maximum number of layers in the model + mmax = nlayer +c is the model flat (nsph = 0) or sphere (nsph = 1) + nsph = iflsph + +c----- +c save current values + do 39 i=1,mmax + b(i) = vsm(i) + a(i) = vpm(i) + d(i) = thkm(i) + rho(i) = rhom(i) +c print *,d(i), b(i) + 39 continue + + if(iwave.eq.1)then + idispl = kmax + idispr = 0 + elseif(iwave.eq.2)then + idispl = 0 + idispr = kmax + endif + + iverb(1) = 0 + iverb(2) = 0 +c ---- constant value + sone0 = 1.500 +c ---- phase velocity increment for searching root + ddc0 = 0.005 +c ---- frequency increment (%) for calculating group vel. using g = dw/dk = dw/d(w/c) + h0 = 0.005 +c ---- period range is:ie for calculation of dispersion + +c----- +c check for water layer +c----- + llw=1 + if(b(1).le.0.0) llw=2 + twopi=2.d0*3.141592653589793d0 + one=1.0d-2 + if(nsph.eq.1) call sphere(0,0,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) + JMN = 1 + betmx=-1.e20 + betmn=1.e20 +c----- +c find the extremal velocities to assist in starting search +c----- + do 20 i=1,mmax + if(b(i).gt.0.01 .and. b(i).lt.betmn)then + betmn = b(i) + jmn = i + jsol = 1 + elseif(b(i).le.0.01 .and. a(i).lt.betmn)then + betmn = a(i) + jmn = i + jsol = 0 + endif + if(b(i).gt.betmx) betmx=b(i) + 20 continue +cc WRITE(6,*)'betmn, betmx:',betmn, betmx +c if(idispl.gt.0)then +cc open(1,file='tmpsrfi.06',form='unformatted', +cc 1 access='sequential') +cc rewind 1 +c read(*,*) lovdispfile +c open(1, file = lovdispfile); +c endif +c if(idispr.gt.0)then +cc open(2,file='tmpsrfi.07',form='unformatted', +cc 1 access='sequential') +cc rewind 2 +c read(*,*) raydispfile +c open(2, file = raydispfile); +c endif + do 2000 ifunc=1,2 + if(ifunc.eq.1.and.idispl.le.0) go to 2000 + if(ifunc.eq.2.and.idispr.le.0) go to 2000 + if(nsph.eq.1) call sphere(ifunc,1,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) + ddc = ddc0 + sone = sone0 + h = h0 +c read(*,*) kmax,mode,ddc,sone,igr,h +c write(*,*) kmax,mode,ddc,sone,igr,h +c read(*,*) (t(i),i=1,kmax) +c write(*,*) (t(i),i=1,kmax) +cc write(ifunc,*) mmax,nsph +cc write(ifunc,*) (btp(i),i=1,mmax) +cc write(ifunc,*) (dtp(i),i=1,mmax) +cc do 420 i=1,mmax +cc write(ifunc,*) d(i),a(i),b(i),rho(i) +cc 420 continue +c write(ifunc,*) kmax,igr,h + if(sone.lt. 0.01) sone=2.0 + onea=dble(sone) +c----- +c get starting value for phase velocity, +c which will correspond to the +c VP/VS ratio +c----- + if(jsol.eq.0)then +c----- +c water layer +c----- + cc1 = betmn + else +c----- +c solid layer solve halfspace period equation +c----- + call gtsolh(a(jmn),b(jmn),cc1) + endif +c----- +c back off a bit to get a starting value at a lower phase velocity +c----- + cc1=.95*cc1 + CC1=.90*CC1 + cc=dble(cc1) + dc=dble(ddc) + dc = dabs(dc) + c1=cc + cm=cc + do 450 i=1,kmax + cb(i)=0.0d0 + c(i)=0.0d0 + 450 continue + ift=999 + do 1800 iq=1,mode + is = 1 + ie = kmax +c read(*,*) is,ie +c write(*,*) 'is =', is, ', ie = ', ie + itst=ifunc + do 1600 k=is,ie + if(k.ge.ift) go to 1700 + t1=dble(t(k)) + if(igr.gt.0)then + t1a=t1/(1.+h) + t1b=t1/(1.-h) + t1=dble(t1a) + else + t1a=sngl(t1) + tlb=0.0 + endif +c----- +c get initial phase velocity estimate to begin search +c +c in the notation here, c() is an array of phase velocities +c c(k-1) is the velocity estimate of the present mode +c at the k-1 period, while c(k) is the phase velocity of the +c previous mode at the k period. Since there must be no mode +c crossing, we make use of these values. The only complexity +c is that the dispersion may be reversed. +c +c The subroutine getsol determines the zero crossing and refines +c the root. +c----- + if(k.eq.is .and. iq.eq.1)then + c1 = cc + clow = cc + ifirst = 1 + elseif(k.eq.is .and. iq.gt.1)then + c1 = c(is) + one*dc + clow = c1 + ifirst = 1 + elseif(k.gt.is .and. iq.gt.1)then + ifirst = 0 +c clow = c(k) + one*dc +c c1 = c(k-1) -onea*dc + clow = c(k) + one*dc + c1 = c(k-1) + if(c1 .lt. clow)c1 = clow + elseif(k.gt.is .and. iq.eq.1)then + ifirst = 0 + c1 = c(k-1) - onea*dc + clow = cm + endif +c----- +c bracket root and refine it +c----- + call getsol(t1,c1,clow,dc,cm,betmx,iret,ifunc,ifirst,d,a,b,rho,rtp,dtp,btp,mmax,llw) + if(iret.eq.-1)goto 1700 + c(k) = c1 +c----- +c for group velocities compute near above solution +c----- + if(igr.gt.0) then + t1=dble(t1b) + ifirst = 0 + clow = cb(k) + one*dc + c1 = c1 -onea*dc + call getsol(t1,c1,clow,dc,cm,betmx,iret,ifunc,ifirst,d,a,b,rho,rtp,dtp,btp,mmax,llw) +c----- +c test if root not found at slightly larger period +c----- + if(iret.eq.-1)then + c1 = c(k) + endif + cb(k)=c1 + else + c1 = 0.0d+00 + endif + cc0 = sngl(c(k)) + cc1 = sngl(c1) + if(igr.eq.0) then +c ----- output only phase velocity +c write(ifunc,*) itst,iq,t(k),cc0,0.0 + cg(k) = cc0 + else +c ----- calculate group velocity and output phase and group velocities + gvel = (1/t1a-1/t1b)/(1/(t1a*cc0)-1/(t1b*cc1)) + cg(k) = gvel +c write(ifunc,*) itst,iq,t(k),(cc0+cc1)/2,gvel +c ----- print *, itst,iq,t(k),t1a,t1b,cc0,cc1,gvel + endif + 1600 continue + go to 1800 + 1700 if(iq.gt.1) go to 1750 + if(iverb(ifunc).eq.0)then + iverb(ifunc) = 1 + write(LOT,*)'improper initial value in disper - no zero found' + write(*,*)'WARNING:improper initial value in disper - no zero found' + write(LOT,*)'in fundamental mode ' + write(LOT,*)'This may be due to low velocity zone ' + write(LOT,*)'causing reverse phase velocity dispersion, ' + write(LOT,*)'and mode jumping.' + write(LOT,*)'due to looking for Love waves in a halfspace' + write(LOT,*)'which is OK if there are Rayleigh data.' + write(LOT,*)'If reverse dispersion is the problem,' + write(LOT,*)'Get present model using OPTION 28, edit sobs.d,' + write(LOT,*)'Rerun with onel large than 2' + write(LOT,*)'which is the default ' +c----- +c if we have higher mode data and the model does not find that +c mode, just indicate (itst=0) that it has not been found, but +c fill out file with dummy results to maintain format - note +c eigenfunctions will not be found for these values. The subroutine +c 'amat' in 'surf' will worry about this in building up the +c input file for 'surfinv' +c----- + write(LOT,*)'ifunc = ',ifunc ,' (1=L, 2=R)' + write(LOT,*)'mode = ',iq-1 + write(LOT,*)'period= ',t(k), ' for k,is,ie=',k,is,ie + write(LOT,*)'cc,cm = ',cc,cm + write(LOT,*)'c1 = ',c1 + write(LOT,*)'d,a,b,rho (d(mmax)=control ignore)' + write(LOT,'(4f15.5)')(d(i),a(i),b(i),rho(i),i=1,mmax) + write(LOT,*)' c(i),i=1,k (NOTE may be part)' + write(LOT,*)(c(i),i=1,k) + endif +c if(k.gt.0)goto 1750 +c go to 2000 + 1750 ift=k + itst=0 + do 1770 i=k,ie + t1a=t(i) +c write(ifunc,*) itst,iq,t1a,0.0,0.0 + cg(i) = 0.0 + 1770 continue + 1800 continue +c close(ifunc,status='keep') + 2000 continue +c close(3,status='keep') + + end + + + + + + + subroutine gtsolh(a,b,c) +c----- +c starting solution +c----- + real*4 kappa, k2, gk2 + c = 0.95*b + do 100 i=1,5 + gamma = b/a + kappa = c/b + k2 = kappa**2 + gk2 = (gamma*kappa)**2 + fac1 = sqrt(1.0 - gk2) + fac2 = sqrt(1.0 - k2) + fr = (2.0 - k2)**2 - 4.0*fac1*fac2 + frp = -4.0*(2.0-k2) *kappa + 1 +4.0*fac2*gamma*gamma*kappa/fac1 + 2 +4.0*fac1*kappa/fac2 + frp = frp/b + c = c - fr/frp + 100 continue + return + end + + subroutine getsol(t1,c1,clow,dc,cm,betmx,iret,ifunc,ifirst,d,a,b,rho,rtp,dtp,btp,mmax,llw) +c----- +c subroutine to bracket dispersion curve +c and then refine it +c----- +c t1 - period +c c1 - initial guess on low side of mode +c clow - lowest possible value for present mode in a +c reversed direction search +c dc - phase velocity search increment +c cm - minimum possible solution +c betmx - maximum shear velocity +c iret - 1 = successful +c - -1= unsuccessful +c ifunc - 1 - Love +c - 2 - Rayleigh +c ifirst - 1 this is first period for a particular mode +c - 0 this is not the first period +c (this is to define period equation sign +c for mode jumping test) +c----- + parameter (NL=200) + real*8 wvno, omega, twopi + real*8 c1, c2, cn, cm, dc, t1, clow + real*8 dltar, del1, del2, del1st, plmn + save del1st + real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL) + integer llw,mmax +c----- +c to avoid problems in mode jumping with reversed dispersion +c we note what the polarity of period equation is for phase +c velocities just beneath the zero crossing at the +c first period computed. +c----- +c bracket solution +c----- + twopi=2.d0*3.141592653589793d0 + omega=twopi/t1 + wvno=omega/c1 + del1 = dltar(wvno,omega,ifunc,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) + if(ifirst.eq.1)del1st = del1 + plmn = dsign(1.0d+00,del1st)*dsign(1.0d+00,del1) + if(ifirst.eq.1)then + idir = +1 + elseif(ifirst.ne.1 .and. plmn.ge.0.0d+00)then + idir = +1 + elseif(ifirst.ne.1 .and. plmn.lt.0.0d+00)then + idir = -1 + endif +c----- +c idir indicates the direction of the search for the +c true phase velocity from the initial estimate. +c Usually phase velocity increases with period and +c we always underestimate, so phase velocity should increase +c (idir = +1). For reversed dispersion, we should look +c downward from the present estimate. However, we never +c go below the floor of clow, when the direction is reversed +c----- + 1000 continue + if(idir.gt.0)then + c2 = c1 + dc + else + c2 = c1 - dc + endif + if(c2.le.clow)then + idir = +1 + c1 = clow + endif + if(c2.le.clow)goto 1000 + omega=twopi/t1 + wvno=omega/c2 + del2 = dltar(wvno,omega,ifunc,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) + if (dsign(1.0d+00,del1).ne.dsign(1.0d+00,del2)) then + go to 1300 + endif + c1=c2 + del1=del2 +c check that c1 is in region of solutions + if(c1.lt.cm) go to 1700 + if(c1.ge.(betmx+dc)) go to 1700 + go to 1000 +c----- +c root bracketed, refine it +c----- + 1300 call nevill(t1,c1,c2,del1,del2,ifunc,cn,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) + c1 = cn + if(c1.gt.(betmx)) go to 1700 + iret = 1 + return + 1700 continue + iret = -1 + return + end +c +c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +c + subroutine sphere(ifunc,iflag,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) +c----- +c Transform spherical earth to flat earth +c +c Schwab, F. A., and L. Knopoff (1972). Fast surface wave and free +c mode computations, in Methods in Computational Physics, +c Volume 11, +c Seismology: Surface Waves and Earth Oscillations, +c B. A. Bolt (ed), +c Academic Press, New York +c +c Love Wave Equations 44, 45 , 41 pp 112-113 +c Rayleigh Wave Equations 102, 108, 109 pp 142, 144 +c +c Revised 28 DEC 2007 to use mid-point, assume linear variation in +c slowness instead of using average velocity for the layer +c Use the Biswas (1972:PAGEOPH 96, 61-74, 1972) density mapping +c +c ifunc I*4 1 - Love Wave +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) + 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 +c common/para/ mmax,llw,twopi + double precision z0,z1,r0,r1,dr,ar,tmp,twopi + save dhalf + ar=6370.0d0 + dr=0.0d0 + r0=ar + d(mmax)=1.0 + if(iflag.eq.0) then + do 5 i=1,mmax + dtp(i)=d(i) + rtp(i)=rho(i) + 5 continue + do 10 i=1,mmax + dr=dr+dble(d(i)) + r1=ar-dr + z0=ar*dlog(ar/r0) + z1=ar*dlog(ar/r1) + d(i)=z1-z0 +c----- +c use layer midpoint +c----- + TMP=(ar+ar)/(r0+r1) + a(i)=a(i)*tmp + b(i)=b(i)*tmp + btp(i)=tmp + r0=r1 + 10 continue + dhalf = d(mmax) + else + d(mmax) = dhalf + do 30 i=1,mmax + if(ifunc.eq.1)then + rho(i)=rtp(i)*btp(i)**(-5) + else if(ifunc.eq.2)then + rho(i)=rtp(i)*btp(i)**(-2.275) + endif + 30 continue + endif + d(mmax)=0.0 + return + end +c +c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +c + subroutine nevill(t,c1,c2,del1,del2,ifunc,cc,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) +c----- +c hybrid method for refining root once it has been bracketted +c between c1 and c2. interval halving is used where other schemes +c would be inefficient. once suitable region is found neville s +c iteration method is used to find root. +c the procedure alternates between the interval halving and neville +c techniques using whichever is most efficient +c----- +c the control integer nev means the following: +c +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) + 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) + integer llw,mmax +c common/modl/ d,a,b,rho,rtp,dtp,btp +c common/para/ mmax,llw,twopi +c----- +c initial guess +c----- + omega = twopi/t + call half(c1,c2,c3,del3,omega,ifunc,d,a,b,rho,rtp,dtp,btp, + & mmax,llw,twopi,a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + nev = 1 + nctrl=1 + 100 continue + nctrl=nctrl+1 + if(nctrl.ge.100) go to 1000 +c----- +c make sure new estimate is inside the previous values. If not +c perform interval halving +c----- + if(c3 .lt. dmin1(c1,c2) .or. c3. gt.dmax1(c1,c2))then + nev = 0 + call half(c1,c2,c3,del3,omega,ifunc,d,a,b,rho,rtp,dtp,btp, + & mmax,llw,twopi,a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + endif + s13 = del1 - del3 + s32 = del3 - del2 +c----- +c define new bounds according to the sign of the period equation +c----- + if(dsign(1.d+00,del3)*dsign(1.d+00,del1) .lt.0.0d+00)then + c2 = c3 + del2 = del3 + else + c1 = c3 + del1 = del3 + endif +c----- +c check for convergence. A relative error criteria is used +c----- + if(dabs(c1-c2).le.1.d-6*c1) go to 1000 +c----- +c if the slopes are not the same between c1, c3 and c3 +c do not use neville iteration +c----- + if(dsign (1.0d+00,s13).ne.dsign (1.0d+00,s32)) nev = 0 +c----- +c if the period equation differs by more than a factor of 10 +c use interval halving to avoid poor behavior of polynomial fit +c----- + ss1=dabs(del1) + s1=0.01*ss1 + ss2=dabs(del2) + s2=0.01*ss2 + if(s1.gt.ss2.or.s2.gt.ss1 .or. nev.eq.0) then + call half(c1,c2,c3,del3,omega,ifunc,d,a,b,rho,rtp,dtp,btp, + & mmax,llw,twopi,a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + nev = 1 + m = 1 + else + if(nev.eq.2)then + x(m+1) = c3 + y(m+1) = del3 + else + x(1) = c1 + y(1) = del1 + x(2) = c2 + y(2) = del2 + m = 1 + endif +c----- +c perform Neville iteration. Note instead of generating y(x) +c we interchange the x and y of formula to solve for x(y) when +c y = 0 +c----- + do 900 kk = 1,m + j = m-kk+1 + denom = y(m+1) - y(j) + if(dabs(denom).lt.1.0d-10*abs(y(m+1)))goto 950 + x(j)=(-y(j)*x(j+1)+y(m+1)*x(j))/denom + 900 continue + c3 = x(1) + wvno = omega/c3 + del3 = dltar(wvno,omega,ifunc,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) +c & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + nev = 2 + m = m + 1 + if(m.gt.10)m = 10 + goto 951 + 950 continue + call half(c1,c2,c3,del3,omega,ifunc,d,a,b,rho,rtp,dtp,btp, + & mmax,llw,twopi,a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + nev = 1 + m = 1 + 951 continue + endif + goto 100 + 1000 continue + cc = c3 + return + end + + subroutine half(c1,c2,c3,del3,omega,ifunc,d,a,b,rho,rtp,dtp,btp, + & mmax,llw,twopi,a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + implicit double precision (a-h,o-z) + parameter(NL=200) + real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL) + c3 = 0.5*(c1 + c2) + wvno=omega/c3 + del3 = dltar(wvno,omega,ifunc,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) +c & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + return + end +c +c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +c + function dltar(wvno,omega,kk,d,a,b,rho,rtp,dtp,btp,mmax,llw,twop) +c & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) +c control the way to P-SV or SH. +c + implicit double precision (a-h,o-z) + parameter(NL=200) + real*4 d(NL),a(NL),b(NL),rho(NL),rtp(NL),dtp(NL),btp(NL) +c + if(kk.eq.1)then +c love wave period equation + dltar = dltar1(wvno,omega,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) + elseif(kk.eq.2)then +c rayleigh wave period equation + dltar = dltar4(wvno,omega,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) +c & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + endif + end +c +c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +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) + 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 +c common/modl/ d,a,b,rho,rtp,dtp,btp +c common/para/ mmax,llw,twopi +c +c Haskell-Thompson love wave formulation from halfspace +c to surface. +c + beta1=dble(b(mmax)) + rho1=dble(rho(mmax)) + xkb=omega/beta1 + wvnop=wvno+xkb + wvnom=dabs(wvno-xkb) + rb=dsqrt(wvnop*wvnom) + e1=rho1*rb + e2=1.d+00/(beta1*beta1) + mmm1 = mmax - 1 + do 600 m=mmm1,llw,-1 + beta1=dble(b(m)) + rho1=dble(rho(m)) + xmu=rho1*beta1*beta1 + xkb=omega/beta1 + wvnop=wvno+xkb + wvnom=dabs(wvno-xkb) + rb=dsqrt(wvnop*wvnom) + q = dble(d(m))*rb + if(wvno.lt.xkb)then + sinq = dsin(q) + y = sinq/rb + z = -rb*sinq + cosq = dcos(q) + elseif(wvno.eq.xkb)then + cosq=1.0d+00 + y=dble(d(m)) + z=0.0d+00 + else + fac = 0.0d+00 + if(q.lt.16)fac = dexp(-2.0d+0*q) + cosq = ( 1.0d+00 + fac ) * 0.5d+00 + sinq = ( 1.0d+00 - fac ) * 0.5d+00 + y = sinq/rb + z = rb*sinq + endif + e10=e1*cosq+e2*xmu*z + e20=e1*y/xmu+e2*cosq + xnor=dabs(e10) + ynor=dabs(e20) + if(ynor.gt.xnor) xnor=ynor + if(xnor.lt.1.d-40) xnor=1.0d+00 + e1=e10/xnor + e2=e20/xnor + 600 continue + dltar1=e1 + return + end +c +c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +c + function dltar4(wvno,omga,d,a,b,rho,rtp,dtp,btp,mmax,llw,twopi) +c & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) +c find P-SV dispersion values. +c + parameter (NL=200,NP=60) + 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) +c common/modl/ d,a,b,rho,rtp,dtp,btp +c common/para/ mmax,llw,twopi +c common/ovrflw/ a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz +c + omega=omga + if(omega.lt.1.0d-4) omega=1.0d-4 + wvno2=wvno*wvno + xka=omega/dble(a(mmax)) + xkb=omega/dble(b(mmax)) + wvnop=wvno+xka + wvnom=dabs(wvno-xka) + ra=dsqrt(wvnop*wvnom) + wvnop=wvno+xkb + wvnom=dabs(wvno-xkb) + rb=dsqrt(wvnop*wvnom) + t = dble(b(mmax))/omega +c----- +c E matrix for the bottom half-space. +c----- + gammk = 2.d+00*t*t + gam = gammk*wvno2 + gamm1 = gam - 1.d+00 + rho1=dble(rho(mmax)) + e(1)=rho1*rho1*(gamm1*gamm1-gam*gammk*ra*rb) + e(2)=-rho1*ra + e(3)=rho1*(gamm1-gammk*ra*rb) + e(4)=rho1*rb + e(5)=wvno2-ra*rb +c----- +c matrix multiplication from bottom layer upward +c----- + mmm1 = mmax-1 + do 500 m = mmm1,llw,-1 + xka = omega/dble(a(m)) + xkb = omega/dble(b(m)) + t = dble(b(m))/omega + gammk = 2.d+00*t*t + gam = gammk*wvno2 + wvnop=wvno+xka + wvnom=dabs(wvno-xka) + ra=dsqrt(wvnop*wvnom) + wvnop=wvno+xkb + wvnom=dabs(wvno-xkb) + rb=dsqrt(wvnop*wvnom) + dpth=dble(d(m)) + rho1=dble(rho(m)) + p=ra*dpth + q=rb*dpth + beta=dble(b(m)) +c----- +c evaluate cosP, cosQ,.... in var. +c evaluate Dunkin's matrix in dnka. +c----- + call var(p,q,ra,rb,wvno,xka,xkb,dpth,w,cosp,exa, + & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + call dnka(ca,wvno2,gam,gammk,rho1, + & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + do 200 i=1,5 + cr=0.0d+00 + do 100 j=1,5 + cr=cr+e(j)*ca(j,i) + 100 continue + ee(i)=cr + 200 continue + call normc(ee,exa) + do 300 i = 1,5 + e(i)=ee(i) + 300 continue + 500 continue + if(llw.ne.1) then +c----- +c include water layer. +c----- + xka = omega/dble(a(1)) + wvnop=wvno+xka + wvnom=dabs(wvno-xka) + ra=dsqrt(wvnop*wvnom) + dpth=dble(d(1)) + rho1=dble(rho(1)) + p = ra*dpth + beta = dble(b(1)) + znul = 1.0d-05 + call var(p,znul,ra,znul,wvno,xka,znul,dpth,w,cosp,exa, + & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) + w0=-rho1*w + dltar4 = cosp*e(1) + w0*e(2) + else + dltar4 = e(1) + endif + return + end +c +c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + subroutine var(p,q,ra,rb,wvno,xka,xkb,dpth,w,cosp,exa,a0,cpcq, + & cpy,cpz,cqw,cqx,xy,xz,wy,wz) +c----- +c find variables cosP, cosQ, sinP, sinQ, etc. +c as well as cross products required for compound matrix +c----- +c To handle the hyperbolic functions correctly for large +c arguments, we use an extended precision procedure, +c keeping in mind that the maximum precision in double +c precision is on the order of 16 decimal places. +c +c So cosp = 0.5 ( exp(+p) + exp(-p)) +c = exp(p) * 0.5 * ( 1.0 + exp(-2p) ) +c becomes +c cosp = 0.5 * (1.0 + exp(-2p) ) with an exponent p +c In performing matrix multiplication, we multiply the modified +c cosp terms and add the exponents. At the last step +c when it is necessary to obtain a true amplitude, +c we then form exp(p). For normalized amplitudes at any depth, +c we carry an exponent for the numerator and the denominator, and +c scale the resulting ratio by exp(NUMexp - DENexp) +c +c The propagator matrices have three basic terms +c +c HSKA cosp cosq +c DUNKIN cosp*cosq 1.0 +c +c When the extended floating point is used, we use the +c largest exponent for each, which is the following: +c +c Let pex = p exponent > 0 for evanescent waves = 0 otherwise +c Let sex = s exponent > 0 for evanescent waves = 0 otherwise +c Let exa = pex + sex +c +c Then the modified matrix elements are as follow: +c +c Haskell: cosp -> 0.5 ( 1 + exp(-2p) ) exponent = pex +c cosq -> 0.5 ( 1 + exp(-2q) ) * exp(q-p) +c exponent = pex +c (this is because we are normalizing all elements in the +c Haskell matrix ) +c Compound: +c cosp * cosq -> normalized cosp * cosq exponent = pex + qex +c 1.0 -> exp(-exa) +c----- + implicit double precision (a-h,o-z) +c common/ovrflw/ a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz + exa=0.0d+00 + a0=1.0d+00 +c----- +c examine P-wave eigenfunctions +c checking whether c> vp c=vp or c < vp +c----- + pex = 0.0d+00 + sex = 0.0d+00 + if(wvno.lt.xka)then + sinp = dsin(p) + w=sinp/ra + x=-ra*sinp + cosp=dcos(p) + elseif(wvno.eq.xka)then + cosp = 1.0d+00 + w = dpth + x = 0.0d+00 + elseif(wvno.gt.xka)then + pex = p + fac = 0.0d+00 + if(p.lt.16)fac = dexp(-2.0d+00*p) + cosp = ( 1.0d+00 + fac) * 0.5d+00 + sinp = ( 1.0d+00 - fac) * 0.5d+00 + w=sinp/ra + x=ra*sinp + endif +c----- +c examine S-wave eigenfunctions +c checking whether c > vs, c = vs, c < vs +c----- + if(wvno.lt.xkb)then + sinq=dsin(q) + y=sinq/rb + z=-rb*sinq + cosq=dcos(q) + elseif(wvno.eq.xkb)then + cosq=1.0d+00 + y=dpth + z=0.0d+00 + elseif(wvno.gt.xkb)then + sex = q + fac = 0.0d+00 + if(q.lt.16)fac = dexp(-2.0d+0*q) + cosq = ( 1.0d+00 + fac ) * 0.5d+00 + sinq = ( 1.0d+00 - fac ) * 0.5d+00 + y = sinq/rb + z = rb*sinq + endif +c----- +c form eigenfunction products for use with compound matrices +c----- + exa = pex + sex + a0=0.0d+00 + if(exa.lt.60.0d+00) a0=dexp(-exa) + cpcq=cosp*cosq + cpy=cosp*y + cpz=cosp*z + cqw=cosq*w + cqx=cosq*x + xy=x*y + xz=x*z + wy=w*y + wz=w*z + qmp = sex - pex + fac = 0.0d+00 + if(qmp.gt.-40.0d+00)fac = dexp(qmp) + cosq = cosq*fac + y=fac*y + z=fac*z + return + end +c +c +c + subroutine normc(ee,ex) +c This routine is an important step to control over- or +c underflow. +c The Haskell or Dunkin vectors are normalized before +c the layer matrix stacking. +c Note that some precision will be lost during normalization. +c + implicit double precision (a-h,o-z) + dimension ee(5) + ex = 0.0d+00 + t1 = 0.0d+00 + do 10 i = 1,5 + if(dabs(ee(i)).gt.t1) t1 = dabs(ee(i)) + 10 continue + if(t1.lt.1.d-40) t1=1.d+00 + do 20 i =1,5 + t2=ee(i) + t2=t2/t1 + ee(i)=t2 + 20 continue +c----- +c store the normalization factor in exponential form. +c----- + ex=dlog(t1) + return + end +c +c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +c + subroutine dnka(ca,wvno2,gam,gammk,rho, + & a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz) +c Dunkin's matrix. +c + implicit double precision (a-h,o-z) + dimension ca(5,5) +c common/ ovrflw / a0,cpcq,cpy,cpz,cqw,cqx,xy,xz,wy,wz + data one,two/1.d+00,2.d+00/ + gamm1 = gam-one + twgm1=gam+gamm1 + gmgmk=gam*gammk + gmgm1=gam*gamm1 + gm1sq=gamm1*gamm1 + rho2=rho*rho + a0pq=a0-cpcq + ca(1,1)=cpcq-two*gmgm1*a0pq-gmgmk*xz-wvno2*gm1sq*wy + ca(1,2)=(wvno2*cpy-cqx)/rho + ca(1,3)=-(twgm1*a0pq+gammk*xz+wvno2*gamm1*wy)/rho + ca(1,4)=(cpz-wvno2*cqw)/rho + ca(1,5)=-(two*wvno2*a0pq+xz+wvno2*wvno2*wy)/rho2 + ca(2,1)=(gmgmk*cpz-gm1sq*cqw)*rho + ca(2,2)=cpcq + ca(2,3)=gammk*cpz-gamm1*cqw + ca(2,4)=-wz + ca(2,5)=ca(1,4) + ca(4,1)=(gm1sq*cpy-gmgmk*cqx)*rho + ca(4,2)=-xy + ca(4,3)=gamm1*cpy-gammk*cqx + ca(4,4)=ca(2,2) + ca(4,5)=ca(1,2) + ca(5,1)=-(two*gmgmk*gm1sq*a0pq+gmgmk*gmgmk*xz+ + * gm1sq*gm1sq*wy)*rho2 + ca(5,2)=ca(4,1) + ca(5,3)=-(gammk*gamm1*twgm1*a0pq+gam*gammk*gammk*xz+ + * gamm1*gm1sq*wy)*rho + ca(5,4)=ca(2,1) + ca(5,5)=ca(1,1) + t=-two*wvno2 + ca(3,1)=t*ca(5,3) + ca(3,2)=t*ca(4,3) + ca(3,3)=a0+two*(cpcq-ca(1,1)) + ca(3,4)=t*ca(2,3) + ca(3,5)=t*ca(1,3) + return + end