C
C Program GPSTEP to calculate the discretization steps C for the decomposition of the wavefield into optimized C 2-D Gaussian packets C C Version: 6.10 C Date: 2007, June 14 C C Coded by: Karel Zacek C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C C....................................................................... C C Description of data files: C C Main input data read from external interactive device (*): C The data consist of character strings, read by list directed (free C format) input. The strings have thus to be enclosed in C apostrophes. The interactive * external unit may be redirected to C the file containing the data. C (1) 'SEP'/ C 'SEP'...Name of the file with input parameters. C Description of file SEP C If blank, default values of the corresponding data are C considered. C Default: 'SEP'=' '. C C C Input file 'SEP' in the SEP format: C The file has the form of a SEP parameter file. C For the description of the SEP format refer to file C 'sep.for'. C Names of input files: C FGBYMIN='string'... Name of the input file containing the C minimum value of the optimized and smoothed C initial parameter of Gaussian beams Y0. C Default: FGBYMIN='fgbymin.out' C FGBYMAX='string'... Name of the input file containing the C maximum value of the optimized and smoothed C initial parameter of Gaussian beams Y0. C Default: FGBYMAX='fgbymax.out' C FGBRMIN='string'... Name of the input file containing the C minimum value of the optimized and smoothed C initial parameter of Gaussian beams R0. C Default: FGBRMIN='fgbrmin.out' C FGBRMAX='string'... Name of the input file containing the C maximum value of the optimized and smoothed C initial parameter of Gaussian beams R0. C Default: FGBRMAX='fgbrmax.out' C MODEL='string'... Input data file describing the model, it is C described in the subroutine file 'model.for'. C Description of file MODEL C Default: 'MODEL'='gpm-mod.dat' C Data specifying the dimensions of the input grid: C NT=positive integer... Number of gridpoints along the time axis T C of the decomposed wavefield. C Default: NT=750 C NX=positive integer... Number of gridpoints along the profile C axis X of the decomposed wavefield. C Default: NX=116 C DT=real... Grid interval along the time axis T of the decomposed C wavefield. C Default: DT=0.004 C DX=real... Grid interval along the profile axis X of C the decomposed wavefield. C Default: DX=25.0 C Name of input parameters: C SHOT=positive integer... Specifies the position of shot. C Default: SHOT=1 C PARK=real... Parameter ksi with values in the interval (0,1). C Default: PARK=0.5 C PARN=real... Parameters ny with values in the interval (0,1). C Default: PARN=0.5 C KAPPA=real... Parameter kappa has the meaning of the ratio of C the step between two subsequent beams to their C characteristic half-width. C Default: KAPPA=1.253314137 C VSURF=real... Velocity determined according to maximum velocity C in the depth of the receivers. C Default: VSURF=2000. C FMAX=real... Maximum frequency of the Gaussian packets. C Default: FMAX=1. C FMIN=real... Minimum frequency of the Gaussian packets. C Default: FMIN=1. C FREF=real... Reference frequencies of the Gaussian packets. C Default: FREF=1. C Name of the output file: C GPPAR='string'... Name of the output file with parameters C K0, K0~, N0 - real part, N0 - imaginary part, C minimum of optimum initial parameter Y0. C Default: GPPAR='gpm-par.out' C GPGRID='string'... Name of the output file containing C discretization steps OXR,OPR,OTR,OOMEGA,DXR,DPR,DTR, C DOMEGA,NXR,NPR,NTR,NOMEGA C Default: GPGRID='gpgrid.out' C RPAR='string'... Name of the output file containing C ray parameters for back-propagation of the wavefield C from receivers. C Default: RPAR='gpm-rp.dat' C REC='string'... String with the name of the output data file C containing the names and coordinates of receivers. C Default: REC='gpm-rec.dat' C SRC='string'... String with the name of the output data file C containing the names and coordinates of sources. C Default: SRC='gpm-src.dat' C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL RARRAY C RARRAY...File 'forms.for'. EXTERNAL RARRAI C RARRAI...File 'forms.for'. EXTERNAL ERROR C ERROR...File 'error.for'. EXTERNAL RSEP1 C RSEP1...File 'sep.for'. EXTERNAL RSEP3T C RSEP3T...File 'sep.for'. EXTERNAL RSEP3I C RSEP3I...File 'sep.for'. EXTERNAL RSEP3R C RSEP3R...File 'sep.for'. EXTERNAL PARM1 C PARM1...File 'parm.for'. EXTERNAL PARM2 C PARM2...File 'parm.for'. C C----------------------------------------------------------------------- C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INCLUDE 'sep.inc' C sep.inc C INCLUDE 'model.inc' C model.inc C INCLUDE 'auxmod.inc' C auxmod.inc C C----------------------------------------------------------------------- C C Auxiliary storage locations: REAL DX,DT,FMAX,FMIN,OMEGAMAX,VSURF,FREF,OMEGAREF C INTEGER LU1,LU2,LU3,LU4,LU5 INTEGER IRAM(MRAM) INTEGER IXR,IPR,ISHOT INTEGER N,NX,NT,NXR,NTR,NPR,NOMEGA REAL DXR,DTR,DPR,DOMEGA,OXR,OPR,OOMEGA REAL XR,PR,KAPPA,PI,VEL,HUH REAL PRMIN,PRMAX,ANIMIN,ANIMAX,ANRMIN,ANRMAX,PARK,PARN REAL AK0,AK0D,THETA REAL COOR(3),RHO COMPLEX AN,ANDD,AN0,AMPL PARAMETER (LU1=21,LU2=22,LU3=23,LU4=24,LU5=25) PARAMETER (LU6=26,LU7=27,LU8=28,LU9=29,LU10=30,LU11=31) PARAMETER (PI=3.141592654) CHARACTER*80 FILSEP,FGBRMIN,FGBRMAX,FGBYMIN,FGBYMAX CHARACTER*80 GPGRID,GPPAR,RPAR,MODEL,REC,SRC EQUIVALENCE (IRAM,RAM) C C....................................................................... C C Reading main input data: FILSEP=' ' WRITE(*,'(A)')'+GPSTEP: Enter filename : ' READ(*,*) FILSEP IF(FILSEP.EQ.' ') THEN C GPSTEP-01 CALL ERROR('GPSTEP-01: No input file specified') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF WRITE(*,'(A)')'+GPSTEP: Reading, calculating... ' C C Reading all the data from the SEP file into the memory: CALL RSEP1(LU1,FILSEP) C C Reading filenames of the input files CALL RSEP3T('FGBYMIN',FGBYMIN,'fgbymin.out') CALL RSEP3T('FGBYMAX',FGBYMAX,'fgbymax.out') CALL RSEP3T('FGBRMIN',FGBRMIN,'fgbrmin.out') CALL RSEP3T('FGBRMAX',FGBRMAX,'fgbrmax.out') CALL RSEP3T('MODEL',MODEL,'gpm-mod.dat') C C Reading filenames of the output files CALL RSEP3T('GPPAR',GPPAR,'gpm-par.out') CALL RSEP3T('GPGRID',GPGRID,'gpgrid.out') CALL RSEP3T('RPAR',RPAR,'gpm-rp.dat') CALL RSEP3T('REC',REC,'gpm-rec.dat') CALL RSEP3T('SRC',SRC,'gpm-src.dat') C CALL RSEP3I('NT',NT,750) CALL RSEP3I('NX',NX,116) CALL RSEP3R('DT',DT,0.004) CALL RSEP3R('DX',DX,25.0) C CALL RSEP3I('SHOT',ISHOT,1) C CALL RSEP3R('PARK',PARK,0.5) CALL RSEP3R('PARN',PARN,0.5) CALL RSEP3R('KAPPA',KAPPA,1.253314137) CALL RSEP3R('VSURF',VSURF,2000.) CALL RSEP3R('FMAX',FMAX,1.) CALL RSEP3R('FMIN',FMIN,1.) CALL RSEP3R('FREF',FREF,1.) C OPEN(LU2,FILE=FGBYMIN,FORM='FORMATTED',STATUS='OLD') OPEN(LU3,FILE=FGBYMAX,FORM='FORMATTED',STATUS='OLD') OPEN(LU4,FILE=FGBRMIN,FORM='FORMATTED',STATUS='OLD') OPEN(LU5,FILE=FGBRMAX,FORM='FORMATTED',STATUS='OLD') OPEN(LU6,FILE=MODEL,FORM='FORMATTED',STATUS='OLD') C OPEN(LU7,FILE=GPPAR,FORM='FORMATTED') OPEN(LU8,FILE=GPGRID,FORM='FORMATTED') OPEN(LU9,FILE=RPAR,FORM='FORMATTED') OPEN(LU10,FILE=REC,FORM='FORMATTED') OPEN(LU11,FILE=SRC,FORM='FORMATTED') C READ(LU2,*) ANIMIN READ(LU3,*) ANIMAX READ(LU4,*) ANRMIN READ(LU5,*) ANRMAX C OMEGAMAX=2.0*PI*FMAX OOMEGA=2.0*PI*FMIN OMEGAREF=2.0*PI*FREF IF(OOMEGA.LT.(1.0/(DT*(NT-1)))) THEN OOMEGA= 1.0/(DT*(NT-1)) END IF C C AK0=PARK*OOMEGA*ANIMIN*VSURF**2 C AK0D=AK0*(1.0-PARN)/PARN C AN0=PARN*2.0*CMPLX(ANRMIN,ANIMIN) AK0=0.5*OOMEGA*ANIMIN*VSURF**2 AK0D=AK0 AN0=CMPLX(ANRMIN,ANIMIN) C WRITE(*,'(A,1(G16.6,1X))') 'K0 =',AK0 WRITE(*,'(A,1(G16.6,1X))') 'K0D =',AK0D WRITE(*,'(A,1(G16.6,1X))') 'N0R =',REAL(AN0) WRITE(*,'(A,1(G16.6,1X))') 'N0I =',AIMAG(AN0) WRITE(*,'(A,1(G16.6,1X))') 'NMIN =',ANIMIN WRITE(LU7,'(1(G16.6,1X))') AK0 WRITE(LU7,'(1(G16.6,1X))') AK0D WRITE(LU7,'(1(G16.6,1X))') REAL(AN0) WRITE(LU7,'(1(G16.6,1X))') AIMAG(AN0) WRITE(LU7,'(1(G16.6,1X))') ANIMIN C IF(OOMEGA.LT.OMEGAMAX) THEN DOMEGA=KAPPA*SQRT(2.0*AK0*AK0D/(AK0+AK0D)) N=NINT((OMEGAMAX-OOMEGA)/DOMEGA) NOMEGA=N+2 DOMEGA=(OMEGAMAX-OOMEGA)/(N+1) ELSE OOMEGA=OMEGAMAX DOMEGA=1.0 NOMEGA=1 END IF WRITE(*,'(A,1(G16.6,1X))') 'DF =',DOMEGA/(2.0*PI) WRITE(*,'(A,I6)') 'NF =',NOMEGA C C DXR=KAPPA*SQRT(-2.0*(1.0-PARN)/OMEGAREF*AIMAG( C * (CMPLX(1.0,0.0)-OOMEGA/OMEGAREF* C * CMPLX(0.0,ANIMIN)/CMPLX(ANRMAX,ANIMAX)*PARK)/ C * CMPLX(ANRMAX,ANIMAX))) DXR=KAPPA*SQRT(-AIMAG(0.5/(OMEGAREF*AN0))) N=NINT((NX-1)*DX/DXR) NXR=N+2 DXR=(NX-1)*DX/(N+1) WRITE(*,'(A,1(G16.6,1X))') 'DXR =',DXR WRITE(*,'(A,I6)') 'NXR =',NXR C OPR=-1.0/VSURF DPR=KAPPA*SQRT(AIMAG(AN0)/OMEGAREF) N=NINT(2.0/(VSURF*DPR)) NPR=N+2 DPR=2.0/(VSURF*(N+1)) WRITE(*,'(A,1(G16.6,1X))') 'DP =',DPR WRITE(*,'(A,I6)') 'NP =',NPR C C DTR=KAPPA*SQRT(-2.0*PARN*AIMAG( C * (CMPLX(1.0,0.0)-CMPLX(0.0,ANIMIN)/CMPLX(ANRMIN,ANIMIN)*PARK) C * /CMPLX(0.0,AK0))) DTR=KAPPA*SQRT(0.5/AK0) N=NINT((NT-1)*DT/DTR) NTR=N+2 DTR=(NT-1)*DT/(N+1) WRITE(*,'(A,1(G16.6,1X))') 'DTR =',DTR WRITE(*,'(A,I6)') 'NTR =',NTR C OXR=150.0+ISHOT*25.0 OTR=0.0 C WRITE(LU8,'(1(G16.6,1X))') OXR WRITE(LU8,'(1(G16.6,1X))') OPR WRITE(LU8,'(1(G16.6,1X))') OTR WRITE(LU8,'(1(G16.6,1X))') OOMEGA WRITE(LU8,'(1(G16.6,1X))') DXR WRITE(LU8,'(1(G16.6,1X))') DPR WRITE(LU8,'(1(G16.6,1X))') DTR WRITE(LU8,'(1(G16.6,1X))') DOMEGA WRITE(LU8,'(I6)') NXR WRITE(LU8,'(I6)') NPR WRITE(LU8,'(I6)') NTR WRITE(LU8,'(I6)') NOMEGA C C CALL PARM1(LU6,1) CALL MODEL1(LU6) WRITE(LU9,'(A)') '''RPAR''' WRITE(LU9,'(A)') '0 0 0 /' WRITE(LU10,'(A)') '/' DO 20 IXR=1,NXR XR=OXR+(IXR-1)*DXR COOR(1)=0.0 COOR(2)=XR COOR(3)=12.0 WRITE(LU10,'(A,3(G16.6,1X),A)') * '''REC'' ', COOR(1),COOR(2),COOR(3),' /' CALL PARM2(1,COOR,UP,US,RHO,QP,QS) VEL=UP(1) C VEL=2000.0 DO 10 IPR=1,NPR HUH=VEL*(OPR+(IPR-1)*DPR) IF(HUH.LT.-1.0) THEN THETA=-PI/2.0 ELSE IF(HUH.GT.1.0) THEN THETA=PI/2.0 ELSE THETA=ASIN(HUH) END IF WRITE(LU9,'(A,1(G16.6,1X),A)')'0 ',THETA,' 0 2 1 0 0.1 0.1 /' 10 CONTINUE WRITE(LU9,'(A)') '/' 20 CONTINUE WRITE(LU10,'(A)') '/' C XR=3000.0+(ISHOT-1)*25.0 COOR(1)=0.0 COOR(2)=XR COOR(3)=12.0 WRITE(LU11,'(A)') '/' WRITE(LU11,'(A,3(G16.6,1X),A)') * '''SRC'' ', COOR(3),COOR(2),COOR(1),' /' WRITE(LU11,'(A)') '/' C CLOSE(LU2) CLOSE(LU3) CLOSE(LU4) CLOSE(LU5) CLOSE(LU6) CLOSE(LU7) CLOSE(LU8) CLOSE(LU9) WRITE(*,'(A)') '+GPSTEP: Done. ' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'length.for' C length.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'model.for' C model.for INCLUDE 'srfc.for' C srfc.for INCLUDE 'metric.for' C metric.for INCLUDE 'parm.for' C parm.for INCLUDE 'val.for' C val.for INCLUDE 'fit.for' C fit.for C C======================================================================= C