C
C Program GRDBORN to calculate the Born approximation of the wavefield C at specified receivers C C Version 6.60 C Date: 2012, June 5 C C Coded by: Libor Sachl C Department of Geophysics, Charles University in Prague, C E-mail: sachl@karel.troja.mff.cuni.cz C C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Names of the input files describing the medium perturbations: C The files have form GRD (GRiDded data). C VPPER='string'... Name of the file with the grid values of the C perturbations of the product of the density and the square C of P-wave velocity. C Default: VPPER=' ' C VSPER='string'... Name of the file with the grid values of the C perturbations of the product of the density and the square C of S-wave velocity. C Default: VSPER=' ' C RHOPER='string'... Name of the file with the grid values of the C density perturbations. C Default: RHOPER=' ' C Name of the input file specifying the receiver Green functions: C BORN='string'... Name of the file containing the names of the C receivers and the names of the corresponding files with C the quantities required for the Born approximation. C The multivalued quantities are discretized on the grid, C are related to the waves arriving to the receiver from the C gridpoints, and are stored in the files of form C MGRD (Multivalued GRiDded data). C Description of file BORN C No default, BORN must be specified and cannot be blank. C IREC=positive integer... Index of the line in file BORN C corresponding to the first receiver. C Default: IREC=1 C NREC=positive integer... Number of receivers read from file BORN. C Receivers IREC to IREC+NREC-1 are thus considered. C Default: NREC=1 C Names of the input files describing the incident wavefield: C The files have form MGRD (Multivalued GRiDded data). C NUM='string'... Name of the file with N1*N2*N3 array of integer C values, containing the number of arrivals at each C gridpoint of the grid. C No default, NUM must be specified and cannot be blank. C MTT='string'... Name of the file with the array of values, C containing for each gridpoint the travel times C of all determined arrivals. C No default, MTT must be specified and cannot be blank. C MP1='string', MP2='string', MP3='string'... Names of the files C with the array of values, containing the first, second and C third component of the slowness vector (i.e. with the C derivatives of travel time with respect to the coordinate C of the point). C No default, MP1, MP2, MP3 must be specified and cannot be C blank. C MTTXX='string'... Name of the file with the grid values of the C second derivative of travel time with respect to the C coordinate perpendicular to the 2D grid. C This file is used for computing the Born approximation of C a 3D wavefield in a 2D grid covering a 2D velocity model, C with the source and the receivers situated in the symmetry C plane of the 2D model. C MTTXX is not used and need not be specified for a 3D grid. C No default, MTTXX must be specified and cannot be blank C for a 2D grid. C AUR1='string', AUI1='string', C AUR2='string', AUI2='string', C AUR3='string', AUI3='string'... Names of the files with the array C of values containing the real or imaginary part of the C amplitude of the wave field. C No default, AUR1, AUI1, AUR2, AUI2, AUR3, AUI3 must be C specified and cannot be blank. C Input files with source and receiver coordinates: C SRC='string'... String with the name of the input data file C with the name and coordinates of the source. C The point names cannot be longer than 6 characters. C Description of file SRC C No default, SRC must be specified and cannot be blank. C ISRC=positive integer... Index of the point in file SRC C corresponding to the source. C Default: ISRC=1 C REC='string'... String with the name of the input data file C with the coordinates corresponding to the receiver names. C The point names cannot be longer than 6 characters. C Description of file REC C No default, REC must be specified and cannot be blank. C Output file: C RF='string'... Name of the output file with the real and imaginary C parts of the response function corresponding to the Born C approximation. C Description of file RF. C Default: RF='rf.out' C Data specifying dimensions of the input grid: C N1=positive integer... Number of gridpoints along the X1 axis. C Default: N1=1 C N2=positive integer... Number of gridpoints along the X2 axis. C Default: N2=1 C N3=positive integer... Number of gridpoints along the X3 axis. C Default: N3=1 C D1=real... Grid spacing along the X1 axis. C Default: D1=1. C D2=real... Grid spacing along the X2 axis. C Default: D2=1. C D3=real... Grid spacing along the X3 axis. C Default: D3=1. C Data describing the frequency domain common with program 'ss.for': C These data are used to generate the optimum default values of C input parameters DF, OF and NF compatible with data for C 'ss.for'. C DT=real... Time step to digitize the source time function and C seismograms. C Default: DT=1. C NFFT=integer... Number of the time samples for the fast Fourier C transform. Will be used to convert the time step to the C frequency step. C NFFT must be an integer power of 2. C Default: NFFT=1 C FMIN=real... The lowest frequency. C Default: FMIN=0 (mostly sufficient) C FMAX=real... The highest frequency. C Default: FMAX=0.5/DT (Nyquist frequency, mostly C sufficient) C Optional data for alternative specification of the frequency domain: C DF=real... Frequency step. C Default: DF=1/(NFFT*DT) (recommended) C OF=real... Born approximation is calculated for NF frequencies C OF,OF+DF,OF+2*DF,...,OF+(NF-1)*DF. C Default: OF=DF*NINT(FMIN/DF) (i.e. FMIN rounded to the C nearest multiple of the frequency step, recommended) C NF=integer... Number of frequencies. C Default: NF=NINT((FMAX-OF)/DF)+1 (recommended) C Hint: Do not specify DF, OF and NF. Use the data for program C 'ss.for' instead and leave the default values of DF, OF and NF C if you are going to generate synthetic seismograms by 'ss.for'. C Data specifying the cosine window used in the Born approximation. C The length of the cosine tapper is L. The cosine window C vanishes at the distance greater than half the grid C interval from each side of the grid. C CSWIN=real... The global specification of parameter named L in the C specification above for cosine window. If none of C parameters CSWIN1,CSWIN2,CSWIN3,CSWIN4,CSWIN5,CSWIN6 is C specified, L is set to this value. C Default: CSWIN=0 C CSWIN1=real... Parameter named L in the specification above for C the cosine window in the negative direction of the 1-st C axis. C Default: CSWIN1=CSWIN C CSWIN2=real... Parameter named L in the specification above for C the cosine window in the positive direction of the 1-st C axis. C Default: CSWIN2=CSWIN C CSWIN3=real... Parameter named L in the specification above for C the cosine window in the negative direction of the 2-nd C axis. C Default: CSWIN3=CSWIN C CSWIN4=real... Parameter named L in the specification above for C the cosine window in the positive direction of the 2-nd C axis. C Default: CSWIN4=CSWIN C CSWIN5=real... Parameter named L in the specification above for C the cosine window in the negative direction of the 3-rd C axis. C Default: CSWIN5=CSWIN C CSWIN6=real... Parameter named L in the specification above for C the cosine window in the positive direction of the 3-rd C axis. C Default: CSWIN6=CSWIN C C C Input file BORN containing names of the receivers and the names of the C files containing the quantities necessary for the Born approximation: C The quantities are discretized on the grid and are related to the C waves arriving to the receiver from each gridpoint of the grid. C The files have form MGRD (Multivalued GRiDded data). C (1) For each receiver data (1.1): C (1.1) TXTREC,NUMG,MTTG,MPG1,MPG2,MPG3,MTTXXG, C AMPR11,AMPI11,AMPR21,AMPI21,AMPR31,AMPI31, C AMPR12,AMPI12,AMPR22,AMPI22,AMPR32,AMPI32, C AMPR13,AMPI13,AMPR23,AMPI23,AMPR33,AMPI33,/ C TXTREC... Name of the receiver. C The name cannot be longer than 6 characters. C NUMG... Name of the file with N1*N2*N3 array of integer C values containing the number of arrivals. C MTTG... Name of the file containing travel times. C MPG1,MPG2,MPG3... 3 names of the files containg the first, C second and third component of the slowness vector C (i.e. with the derivatives of travel time with respect to C the coordinates of the point). C MTTXXG... Name of the file with the grid values of the second C derivative of travel time with respect to the coordinate C perpendicular to the 2D grid. C This file is used for computing the Born approximation of C a 3D wavefield in a 2D grid covering a 2D velocity model, C with the source and the receivers situated in the symmetry C plane of the 2D model. C MTTXXG is not used and may be blank for a 3D grid. C MTTXXG cannot be blank for a 2D grid. C AMPR11,AMPI11,AMPR21,AMPI21,AMPR31,AMPI31, C AMPR12,AMPI12,AMPR22,AMPI22,AMPR32,AMPI32, C AMPR13,AMPI13,AMPR23,AMPI23,AMPR33,AMPI33... 18 names of the C files containing the real or imaginary part of the C tensorial amplitude of the Green function. C The filenames have no defaults. They must be specified. C Input file BORN may contain the above data only, without any C additional line. C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C C....................................................................... C C Filenames and parameters: CHARACTER*80 FSEP,FNAME,FILREC,FILSRC CHARACTER*80 NUMG,MTTG,MPG1,MPG2,MPG3,MTTXXG CHARACTER*80 AMPR11,AMPI11,AMPR21,AMPI21,AMPR31,AMPI31, & AMPR12,AMPI12,AMPR22,AMPI22,AMPR32,AMPI32, & AMPR13,AMPI13,AMPR23,AMPI23,AMPR33,AMPI33 INTEGER LU1,LU2,LU3,LU4 PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4) REAL PI,P2I,GIANT PARAMETER (PI=3.14159265,P2I=6.28318530,GIANT=2000000000.) C C Input data: INTEGER N1,N2,N3 INTEGER IREC,NREC INTEGER NPTS,NF REAL D1,D2,D3 REAL TD,FMIN,FMAX,FSTEP C C Data needed for output file INTEGER ISRC CHARACTER*6 TXTREC REAL XS1,XS2,XS3 REAL XR1,XR2,XR3 REAL TTMIN,TTMAX C C Memory allocation: INTEGER INUMU,ILA,IMU,IRHO INTEGER IMTTU,IPU1,IPU2,IPU3,ITTXX,IUR1,IUI1,IUR2,IUI2,IUR3,IUI3 INTEGER IBORN INTEGER INUMG,IMTTG,IPG1,IPG2,IPG3,ITTXXG INTEGER IGR11,IGI11,IGR21,IGI21,IGR31,IGI31,IGR12,IGI12,IGR22 INTEGER IGI22,IGR32,IGI32,IGR13,IGI13,IGR23,IGI23,IGR33,IGI33 C C Auxiliary names for C perturbations of elastic moduli and density: REAL RLA,RMU,RHO C travel times and slowness vectors of the incident wavefield: REAL TTU,PU1,PU2,PU3 C second derivative of travel times of the incident wavefield in C the direction perpendicular to the symmetry plane of the C 2D model. REAL TTXX C real and imaginary parts of amplitudes of incident wavefield: REAL UR1,UR2,UR3 REAL UI1,UI2,UI3 C travel times and slowness vectors of the Green function: REAL TTG,PG1,PG2,PG3 C second derivative of travel times of the Green function in the C direction perpendicular to the symmetry plane of the 2D model. REAL TTXXG C real and imaginary parts of amplitudes of the Green function: REAL GR11,GR21,GR31,GR12,GR22,GR32,GR13,GR23,GR33 REAL GI11,GI21,GI31,GI12,GI22,GI32,GI13,GI23,GI33 C travel time of the incident wavefield plus travel time of C the Green function REAL TT C C Other variables: LOGICAL LINPUT,LRLA,LRMU,LRHO,LVPPER,LVSPER,LLRLA,LLRMU,LLRHO,L2D CHARACTER*80 MSSG,NAMPTS,MRECN CHARACTER*1 TEXT INTEGER I,II,J,JJ,K,IR,JF,ICOUNT INTEGER N12,N23,N123,N,NN,N6F,NWIN1,NWIN2,NWIN3,NWIN4,NWIN5,NWIN6 INTEGER J1,J2,J3,J4,J5,J6,J6F INTEGER I01,I02,I03,I04,I05,I06,I07,I08,I09,I10,I11,I12 REAL D123,ARG,ARGINC,RL,RL1,RL2,RL3,RL4,RL5,RL6,WINDOW REAL AUAG1R,AUAG1I,AUAG2R,AUAG2I,AUAG3R,AUAG3I REAL PGAG1R,PGAG1I,PGAG2R,PGAG2I,PGAG3R,PGAG3I REAL PUAUR,PUAUI REAL PUAG1R,PUAG1I,PUAG2R,PUAG2I,PUAG3R,PUAG3I REAL PGAUR,PGAUI,PUPG REAL DER,DEI,B1R,B1I,B2R,B2I,B3R,B3I,BR,BI REAL F,DOM2 REAL AMAX C C Functions INTEGER LENGTH C C----------------------------------------------------------------------- C C Reading input data C C Reading input SEP parameter file: WRITE(*,'(A)') '+GRDBORN: Enter input filename: ' FSEP=' ' READ(*,*) FSEP IF(FSEP.EQ.' ') THEN C GRDBORN-01 CALL ERROR('GRDBORN-01: SEP file not given.') 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. END IF CALL RSEP1(LU1,FSEP) WRITE(*,'(A)') '+GRDBORN: Working ... ' C C Reading input data specifying dimensions of the input grid: CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) N12=N1*N2 N23=N2*N3 N123=N1*N2*N3 C CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) L2D=.FALSE. IF(N2.EQ.1.AND.N3.EQ.1) THEN C D123=ABS(D1) C WRITE(*,'(A)') '+GRDBORN: 1D Born approximation on a 1D grid.' C GRDBORN-14 CALL ERROR('GRDBORN-14: 1D Born approximation is not coded.') C The Born approximation of 1D, 2D or 3D wavefield on a 1D grid C has not been coded yet. ELSE IF(N1.EQ.1.AND.N3.EQ.1) THEN C D123=ABS(D2) C WRITE(*,'(A)') '+GRDBORN: 1D Born approximation on a 1D grid.' C GRDBORN-15 CALL ERROR('GRDBORN-15: 1D Born approximation is not coded.') C The Born approximation of 1D, 2D or 3D wavefield on a 1D grid C has not been coded yet. ELSE IF(N1.EQ.1.AND.N2.EQ.1) THEN C D123=ABS(D2) C WRITE(*,'(A)') '+GRDBORN: 1D Born approximation on a 1D grid.' C GRDBORN-16 CALL ERROR('GRDBORN-16: 1D Born approximation is not coded.') C The Born approximation of 1D, 2D or 3D wavefield on a 1D grid C has not been coded yet. ELSE IF(N1.EQ.1) THEN D123=ABS(D2*D3) WRITE(*,'(A)') '+GRDBORN: 3D Born approximation on a 2D grid.' L2D=.TRUE. C Note that the Born approximation of 2D wavefield on a 2D grid C has not been coded yet. ELSE IF(N2.EQ.1) THEN D123=ABS(D1*D3) WRITE(*,'(A)') '+GRDBORN: 3D Born approximation on a 2D grid.' L2D=.TRUE. C Note that the Born approximation of 2D wavefield on a 2D grid C has not been coded yet. ELSE IF(N3.EQ.1) THEN D123=ABS(D1*D2) WRITE(*,'(A)') '+GRDBORN: 3D Born approximation on a 2D grid.' L2D=.TRUE. C Note that the Born approximation of 2D wavefield on a 2D grid C has not been coded yet. ELSE D123=ABS(D1*D2*D3) WRITE(*,'(A)') '+GRDBORN: 3D Born approximation on a 3D grid.' END IF C C Checking the size of RAM: C Is array RAM(MRAM) allocated in include file 'ram.inc' big enough C to store content of file specified by NUM and to read medium C perturbations? C You may wish to increase the dimension MRAM in file C ram.inc. CALL CHSP(MRAM,4*N123+1,'perturbations of density') C C Reading data specifying dimensions of the output grid: CALL RSEP3I('IREC',IREC,1) CALL RSEP3I('NREC',NREC,1) C C Reading data specifying the cosine windows. CALL RSEP3R('CSWIN',RL,0.) CALL RSEP3R('CSWIN1',RL1,RL) CALL RSEP3R('CSWIN2',RL2,RL) CALL RSEP3R('CSWIN3',RL3,RL) CALL RSEP3R('CSWIN4',RL4,RL) CALL RSEP3R('CSWIN5',RL5,RL) CALL RSEP3R('CSWIN6',RL6,RL) C C Reading data describing the frequency domain: CALL RSEP3I('NFFT ',NPTS ,1) CALL RSEP3R('DT ',TD ,1.) CALL RSEP3R('FMIN ',FMIN ,0.) CALL RSEP3R('FMAX ',FMAX ,0.5/TD) CALL RSEP3R('DF ',FSTEP ,1./(FLOAT(NPTS)*TD)) CALL RSEP3R('OF ',FMIN ,FSTEP*NINT(FMIN/FSTEP)) CALL RSEP3I('NF ',NF ,NINT((FMAX-FMIN)/FSTEP)+1) C Auxiliary variable N6F=6*NF C C Positions of the first elements of appropriate quantities C in array RAM minus 1. C Quantities:number of arrivals for incident wave field. INUMU=1 IRAM(INUMU)=0 C C Reading the numbers of arrivals for incident wave field CALL CHLI('NUM',' ',LU1,N123,.TRUE.,IRAM(INUMU+1),LINPUT) C C Sequential number of the lastly read arrival at the I-th point: DO 10 I=INUMU+1,INUMU+N123 IRAM(I)=IRAM(I)+IRAM(I-1) 10 CONTINUE C N=IRAM(INUMU+N123) C Reading input data specifying model and computation of positions C of the first elements of appropriate quantities in array RAM C minus 1. C Quantities:perturbations of rho*(v_s)^2, C perturbations of rho*(v_p)^2,density perturbations, C travel times,slowness vectors and amplitudes C of the incident wavefield,Born approximation, C number of arrivals and travel times of the Green C function IMU=INUMU+N123 C CALL CHLR('VSPER',' ',LU1,N123,.FALSE.,RAM(IMU+1),LVSPER) LRMU=LVSPER IF(LRMU) THEN ILA=IMU+N123 ELSE ILA=IMU END IF C CALL CHLR('VPPER',' ',LU1,N123,.FALSE.,RAM(ILA+1),LVPPER) IF(LVPPER.OR.LVSPER) THEN LRLA=.TRUE. ELSE LRLA=.FALSE. END IF IF(LRLA) THEN IRHO=ILA+N123 ELSE IRHO=ILA END IF C CALL CHLR('RHOPER',' ',LU1,N123,.FALSE.,RAM(IRHO+1),LRHO) IF(LRHO) THEN IMTTU=IRHO+N123 ELSE IMTTU=IRHO END IF C IPU1=IMTTU+N IPU2=IPU1+N IPU3=IPU2+N C IF(L2D) THEN ITTXX=IPU3+N ELSE ITTXX=IPU3 END IF C IUR1=ITTXX+N IUI1=IUR1+N IUR2=IUI1+N IUI2=IUR2+N IUR3=IUI2+N IUI3=IUR3+N C IBORN=IUI3+N C INUMG=IBORN+N6F+1 IRAM(INUMG)=0 C IMTTG=INUMG+N123 C C Checking the size of RAM: C Is array RAM(MRAM) allocated in include file 'ram.inc' big enough C to contain both input and output grids (without the Green C function)? C You may wish to increase the dimension MRAM in file C ram.inc. CALL CHSP(MRAM,INUMG+N123,'NUMG') C C Reading input data specifying incident wavefield. CALL CHLR('MTT',' ',LU1,N,.TRUE.,RAM(IMTTU+1),LINPUT) C CALL CHLR('MP1',' ',LU1,N,.TRUE.,RAM(IPU1+1),LINPUT) CALL CHLR('MP2',' ',LU1,N,.TRUE.,RAM(IPU2+1),LINPUT) CALL CHLR('MP3',' ',LU1,N,.TRUE.,RAM(IPU3+1),LINPUT) C IF(L2D) THEN CALL CHLR('MTTXX',' ',LU1,N,.TRUE.,RAM(ITTXX+1),LINPUT) DO 17 I=1,N IF(RAM(ITTXX+I).LT.0.) THEN C GRDBORN-12 CALL ERROR('GRDBORN-12: Negative second derivative of ' & // 'travel time in the direction perpendicular ' & // 'to the symmetry plane.') END IF 17 CONTINUE END IF C CALL CHLR('AUR1',' ',LU1,N,.TRUE.,RAM(IUR1+1),LINPUT) CALL CHLR('AUI1',' ',LU1,N,.TRUE.,RAM(IUI1+1),LINPUT) CALL CHLR('AUR2',' ',LU1,N,.TRUE.,RAM(IUR2+1),LINPUT) CALL CHLR('AUI2',' ',LU1,N,.TRUE.,RAM(IUI2+1),LINPUT) CALL CHLR('AUR3',' ',LU1,N,.TRUE.,RAM(IUR3+1),LINPUT) CALL CHLR('AUI3',' ',LU1,N,.TRUE.,RAM(IUI3+1),LINPUT) C C....................................................................... C C Adjustment. Computation of perturbations of elastic moduli lambda C from perturbations of rho*(v_p)^2 and rho*(v_s)^2 : C IF(LVPPER.AND.LVSPER) THEN DO 11 I=1,N123 RAM(ILA+I)=RAM(ILA+I)-2.*RAM(IMU+I) 11 CONTINUE ELSE IF(LVSPER) THEN DO 16 I=1,N123 RAM(ILA+I)=-2.*RAM(IMU+I) 16 CONTINUE END IF C C....................................................................... C C Application of the cosine window on the model boundaries. C The affected gridpoints form a set of gridpoint planes C (two sets of gridpoint planes per an axis). C Number of gridpoint planes in each set: IF(N1.NE.1) THEN NWIN1=INT((RL1+D1/2)/ABS(D1)) NWIN2=INT((RL2+D1/2)/ABS(D1)) ELSE NWIN1=0 NWIN2=0 END IF IF(N2.NE.1) THEN NWIN3=INT((RL3+D2/2)/ABS(D2)) NWIN4=INT((RL4+D2/2)/ABS(D2)) ELSE NWIN3=0 NWIN4=0 END IF IF(N3.NE.1) THEN NWIN5=INT((RL5+D3/2)/ABS(D3)) NWIN6=INT((RL6+D3/2)/ABS(D3)) ELSE NWIN5=0 NWIN6=0 END IF C C 1-st axis C Variables II,J3 and JJ,J4 are pointers belonging to the first or C the second set of gridpoint planes. C Negative direction C Setting initial values of cosine argument, argument increment, C and pointer II. IF (NWIN1.GE.1) THEN ARG=PI*(D1/2.)/RL1 ARGINC=2.*ARG II=1 END IF C Loop over gridpoint values of the cosine window. DO 61 I=1,NWIN1 WINDOW=(1.-COS(ARG))/2. C Pointer J3 is set to the value of pointer II. In the next loop C we will use and alter J3. J3=II C Loop over points with the same value of the cosine window DO 60 J=1,N23 IF(LRHO) THEN J1=IRHO+J3 RAM(J1)=RAM(J1)*WINDOW END IF IF(LRLA) THEN J1=ILA+J3 RAM(J1)=RAM(J1)*WINDOW END IF IF(LRMU) THEN J1=IMU+J3 RAM(J1)=RAM(J1)*WINDOW END IF C Pointer modification. J3=J3+N1 60 CONTINUE C Argument increment ARG=ARG+ARGINC C Pointer modification II=II+1 61 CONTINUE C Positive direction C Setting initial values of cosine argument, argument increment C and pointer JJ. IF (NWIN2.GE.1) THEN ARG=PI*(D1/2.)/RL2 ARGINC=2.*ARG JJ=N1 END IF C Loop over gridpoint values of the cosine window. DO 63 I=1,NWIN2 WINDOW=(1.-COS(ARG))/2. C Pointer J4 is set to the value of pointer JJ. In the next loop C we will use and alter J4. J4=JJ C Loop over points with the same value of the cosine window DO 62 J=1,N23 IF(LRHO) THEN J2=IRHO+J4 RAM(J2)=RAM(J2)*WINDOW END IF IF(LRLA) THEN J2=ILA+J4 RAM(J2)=RAM(J2)*WINDOW END IF IF(LRMU) THEN J2=IMU+J4 RAM(J2)=RAM(J2)*WINDOW END IF C Pointer modification. J4=J4+N1 62 CONTINUE C Argument increment ARG=ARG+ARGINC C Pointer modification JJ=JJ-1 63 CONTINUE C The structure of the code in the following two cases is the same C (3-rd axis) or very similar (2-nd axis), therefore commentary is C briefer C C 2-nd axis C Variables II,J5,J3 and JJ,J6,J4 are pointers belonging to C the first or the second set of gridpoint planes. C Negative direction IF (NWIN3.GE.1) THEN ARG=PI*(D2/2.)/RL3 ARGINC=2.*ARG II=1 END IF C Loop over gridpoint values of the cosine window. DO 66 I=1,NWIN3 WINDOW=(1.-COS(ARG))/2. J5=II C Two loops over points with the same value of the cosine window. DO 65 K=1,N3 J3=J5 DO 64 J=1,N1 IF(LRHO) THEN J1=IRHO+J3 RAM(J1)=RAM(J1)*WINDOW END IF IF(LRLA) THEN J1=ILA+J3 RAM(J1)=RAM(J1)*WINDOW END IF IF(LRMU) THEN J1=IMU+J3 RAM(J1)=RAM(J1)*WINDOW END IF J3=J3+1 64 CONTINUE J5=J5+N12 65 CONTINUE ARG=ARG+ARGINC II=II+N1 66 CONTINUE C Positive direction IF (NWIN4.GE.1) THEN ARG=PI*(D2/2.)/RL4 ARGINC=2.*ARG JJ=N12 END IF C Loop over gridpoint values of the cosine window. DO 69 I=1,NWIN4 WINDOW=(1.-COS(ARG))/2. J6=JJ C Two loops over points with the same value of the cosine window. DO 68 K=1,N3 J4=J6 DO 67 J=1,N1 IF(LRHO) THEN J2=IRHO+J4 RAM(J2)=RAM(J2)*WINDOW END IF IF(LRLA) THEN J2=ILA+J4 RAM(J2)=RAM(J2)*WINDOW END IF IF(LRMU) THEN J2=IMU+J4 RAM(J2)=RAM(J2)*WINDOW END IF J4=J4-1 67 CONTINUE J6=J6+N12 68 CONTINUE ARG=ARG+ARGINC JJ=JJ-N1 69 CONTINUE C C 3-rd axis C Variables II,J3 and JJ,J4 are pointers belonging to the first or C the second set of gridpoint planes. C Negative direction IF (NWIN5.GE.1) THEN ARG=PI*(D3/2.)/RL5 ARGINC=2.*ARG II=1 END IF C Loop over gridpoint values of the cosine window. DO 71 I=1,NWIN5 WINDOW=(1.-COS(ARG))/2. J3=II C Loop over points with the same value of the cosine window DO 70 J=1,N12 IF(LRHO) THEN J1=IRHO+J3 RAM(J1)=RAM(J1)*WINDOW END IF IF(LRLA) THEN J1=ILA+J3 RAM(J1)=RAM(J1)*WINDOW END IF IF(LRMU) THEN J1=IMU+J3 RAM(J1)=RAM(J1)*WINDOW END IF J3=J3+1 70 CONTINUE ARG=ARG+ARGINC II=II+N12 71 CONTINUE C Positive direction IF (NWIN6.GE.1) THEN ARG=PI*(D3/2.)/RL6 ARGINC=2.*ARG JJ=N123 END IF C Loop over gridpoint values of the cosine window. DO 73 I=1,NWIN6 WINDOW=(1.-COS(ARG))/2. J4=JJ C Loop over points with the same value of the cosine window DO 72 J=1,N12 IF(LRHO) THEN J2=IRHO+J4 RAM(J2)=RAM(J2)*WINDOW END IF IF(LRLA) THEN J2=ILA+J4 RAM(J2)=RAM(J2)*WINDOW END IF IF(LRMU) THEN J2=IMU+J4 RAM(J2)=RAM(J2)*WINDOW END IF J4=J4-1 72 CONTINUE ARG=ARG+ARGINC JJ=JJ-N12 73 CONTINUE C C....................................................................... C C Opening input file with the Green function: C CALL RSEP3T('BORN',FNAME ,' ') C C Checking IF(FNAME.EQ.' ') THEN C GRDBORN-02 CALL ERROR('GRDBORN-02: Input file BORN not given.') END IF C OPEN(LU1,FILE=FNAME,STATUS='OLD') DO 12 I=1,IREC-1 READ(LU1,*) TXTREC, & NUMG,MTTG,MPG1,MPG2,MPG3,MTTXXG, & AMPR11,AMPI11,AMPR21,AMPI21,AMPR31,AMPI31, & AMPR12,AMPI12,AMPR22,AMPI22,AMPR32,AMPI32, & AMPR13,AMPI13,AMPR23,AMPI23,AMPR33,AMPI33 12 CONTINUE C C....................................................................... C C Opening the output file and writing parts (1),(2),(3) C C Reading the name of the file with coordinates of the source CALL RSEP3T('SRC',FILSRC,' ') IF(FILSRC.EQ.' ') THEN C GRDBORN-03 CALL ERROR('GRDBORN-03: File with coordinates of source was ' & // 'not specified.') END IF C Reading the position of desired source. CALL RSEP3I('ISRC',ISRC,1) C Reading the coordinates of the source OPEN(LU4,FILE=FILSRC,STATUS='OLD') READ(LU4,*) (TEXT,I=1,20) ICOUNT=0 13 CONTINUE ICOUNT=ICOUNT+1 NAMPTS='$$$$$$$' XS1=0. XS2=0. XS3=0. READ(LU4,*,END=14) NAMPTS,XS1,XS2,XS3 IF(NAMPTS(1:7).EQ.'$$$$$$$') THEN GO TO 14 END IF IF(NAMPTS(7:).NE.' ') THEN C GRDBORN-04 CALL ERROR('GRDBORN-04: Source name exceeds 6 characters.') C Names of points in file given by input parameter PTS cannot C be longer than 6 characters. This limitation is imposed by C the GSE standard. END IF IF(ICOUNT.EQ.ISRC) THEN GO TO 15 END IF GO TO 13 14 CONTINUE C GRDBORN-05 CALL ERROR('GRDBORN-05: Source name and coordinates were not' & // 'specified.') 15 CONTINUE CLOSE(LU4) C Reading the name of the output file and opening that file CALL RSEP3T('RF',FNAME,'rf.out') OPEN(LU3,FILE=FNAME) C Writing. WRITE(LU3,'(A)') '/' WRITE(LU3,'(3(F9.3,1X),A)') XS1,XS2,XS3,'/' WRITE(LU3,'(2(E11.5,1X),I8,1X,A)') FMIN,FSTEP,NF,'/' C C....................................................................... C C Computation C C Reading the name of the file with coordinates of receivers CALL RSEP3T('REC',FILREC,' ') IF(FILREC.EQ.' ') THEN C GRDBORN-06 CALL ERROR('GRDBORN-06: File with coordinates of receivers was ' & // 'not specified.') END IF C Opening the file with coordinates of receivers OPEN(LU4,FILE=FILREC,STATUS='OLD') READ(LU4,*) (TEXT,I=1,20) C C Loop over receivers. DO 50 IR=1,NREC TTMIN=GIANT TTMAX=-GIANT C First of all,elements of array RAM which will contain Born C approximation have to be nullified. DO 20 I=IBORN+1,IBORN+N6F RAM(I)=0. 20 CONTINUE C TXTREC=' ' READ(LU1,*,END=21) TXTREC, & NUMG,MTTG,MPG1,MPG2,MPG3,MTTXXG, & AMPR11,AMPI11,AMPR21,AMPI21,AMPR31,AMPI31, & AMPR12,AMPI12,AMPR22,AMPI22,AMPR32,AMPI32, & AMPR13,AMPI13,AMPR23,AMPI23,AMPR33,AMPI33 C C Checking 21 CONTINUE IF(TXTREC.EQ.' ') THEN WRITE(MRECN,*) IREC+IR-1 C GRDBORN-07 WRITE(MSSG,'(3A)') 'GRDBORN-07: Name of the receiver number ', & MRECN(1:LENGTH(MRECN)), & ' was not specified.' CALL ERROR(MSSG) C Note: NUMG.EQ.' ' is caught by subroutine CHLI2 further on. END IF C C Reading the coordinates of the receiver 22 CONTINUE NAMPTS='$$$$$$$' XR1=0. XR2=0. XR3=0. READ(LU4,*,END=23) NAMPTS,XR1,XR2,XR3 IF(NAMPTS(1:7).EQ.'$$$$$$$') THEN GO TO 23 END IF IF(NAMPTS(7:).NE.' ') THEN C GRDBORN-08 CALL ERROR('GRDBORN-08: Receiver name exceeds ' & // '6 characters.') C Names of points in file given by input parameter PTS cannot C be longer than 6 characters. This limitation is imposed by C the GSE standard. END IF IF(TXTREC.EQ.NAMPTS) THEN GO TO 24 END IF GO TO 22 23 CONTINUE C GRDBORN-09 CALL ERROR('GRDBORN-09: No point name matching the input ' & // 'receiver name.') 24 CONTINUE REWIND(UNIT=LU4) C C Reading the data specifying the number of arrivals for C the Green function. CALL CHLI2('NUMG',NUMG,LU2,N123,.TRUE.,IRAM(INUMG+1),LINPUT) C C Sequential number of the lastly read arrival at the J-th C point: DO 25 J=INUMG+1,INUMG+N123 IRAM(J)=IRAM(J)+IRAM(J-1) 25 CONTINUE C NN=IRAM(INUMG+N123) C C Positions of the first elements of appropriate quantities C in array RAM minus 1. C Quantities:slowness vectors and amplitudes of the Green function C IPG1=IMTTG+NN IPG2=IPG1+NN IPG3=IPG2+NN C IF(L2D) THEN ITTXXG=IPG3+NN ELSE ITTXXG=IPG3 END IF C IGR11=ITTXXG+NN IGI11=IGR11+NN IGR21=IGI11+NN IGI21=IGR21+NN IGR31=IGI21+NN IGI31=IGR31+NN IGR12=IGI31+NN IGI12=IGR12+NN IGR22=IGI12+NN IGI22=IGR22+NN IGR32=IGI22+NN IGI32=IGR32+NN IGR13=IGI32+NN IGI13=IGR13+NN IGR23=IGI13+NN IGI23=IGR23+NN IGR33=IGI23+NN IGI33=IGR33+NN C C Checking CALL CHSP(MRAM,IGI33+NN,'Im(G_33)') C C Reading the Green function CALL CHLR2('MTTG',MTTG,LU2,NN,.TRUE.,RAM(IMTTG+1),LINPUT) C CALL CHLR2('MPG1',MPG1,LU2,NN,.TRUE.,RAM(IPG1+1),LINPUT) CALL CHLR2('MPG2',MPG2,LU2,NN,.TRUE.,RAM(IPG2+1),LINPUT) CALL CHLR2('MPG3',MPG3,LU2,NN,.TRUE.,RAM(IPG3+1),LINPUT) C IF(L2D) THEN CALL CHLR2('MTTXXG',MTTXXG,LU2,NN,.TRUE.,RAM(ITTXXG+1),LINPUT) DO 26 I=1,NN IF(RAM(ITTXXG+I).LT.0.) THEN C GRDBORN-13 CALL ERROR('GRDBORN-13: Negative second derivative of ' & // 'travel time in the direction perpendicular ' & // 'to the symmetry plane.') END IF 26 CONTINUE END IF C CALL CHLR2('AMPR11',AMPR11,LU2,NN,.TRUE.,RAM(IGR11+1),LINPUT) CALL CHLR2('AMPI11',AMPI11,LU2,NN,.TRUE.,RAM(IGI11+1),LINPUT) C CALL CHLR2('AMPR21',AMPR21,LU2,NN,.TRUE.,RAM(IGR21+1),LINPUT) CALL CHLR2('AMPI21',AMPI21,LU2,NN,.TRUE.,RAM(IGI21+1),LINPUT) C CALL CHLR2('AMPR31',AMPR31,LU2,NN,.TRUE.,RAM(IGR31+1),LINPUT) CALL CHLR2('AMPI31',AMPI31,LU2,NN,.TRUE.,RAM(IGI31+1),LINPUT) C CALL CHLR2('AMPR12',AMPR12,LU2,NN,.TRUE.,RAM(IGR12+1),LINPUT) CALL CHLR2('AMPI12',AMPI12,LU2,NN,.TRUE.,RAM(IGI12+1),LINPUT) C CALL CHLR2('AMPR22',AMPR22,LU2,NN,.TRUE.,RAM(IGR22+1),LINPUT) CALL CHLR2('AMPI22',AMPI22,LU2,NN,.TRUE.,RAM(IGI22+1),LINPUT) C CALL CHLR2('AMPR32',AMPR32,LU2,NN,.TRUE.,RAM(IGR32+1),LINPUT) CALL CHLR2('AMPI32',AMPI32,LU2,NN,.TRUE.,RAM(IGI32+1),LINPUT) C CALL CHLR2('AMPR13',AMPR13,LU2,NN,.TRUE.,RAM(IGR13+1),LINPUT) CALL CHLR2('AMPI13',AMPI13,LU2,NN,.TRUE.,RAM(IGI13+1),LINPUT) C CALL CHLR2('AMPR23',AMPR23,LU2,NN,.TRUE.,RAM(IGR23+1),LINPUT) CALL CHLR2('AMPI23',AMPI23,LU2,NN,.TRUE.,RAM(IGI23+1),LINPUT) C CALL CHLR2('AMPR33',AMPR33,LU2,NN,.TRUE.,RAM(IGR33+1),LINPUT) CALL CHLR2('AMPI33',AMPI33,LU2,NN,.TRUE.,RAM(IGI33+1),LINPUT) C WRITE(*,'(A)') '+GRDBORN: Working. ' C Loop over gridpoints. DO 30 I=1,N123 C Auxiliary names and setting local logical variables, which C determine whether perturbations of medium parameters are zero C or not in the current gridpoint. C IF(LRLA) THEN RLA=RAM(ILA+I) IF(RLA.NE.0.) THEN LLRLA=.TRUE. ELSE LLRLA=.FALSE. END IF ELSE LLRLA=.FALSE. END IF C IF(LRMU) THEN RMU=RAM(IMU+I) IF(RMU.NE.0.) THEN LLRMU=.TRUE. ELSE LLRMU=.FALSE. END IF ELSE LLRMU=.FALSE. END IF C IF(LRHO) THEN RHO=RAM(IRHO+I) IF(RHO.NE.0.) THEN LLRHO=.TRUE. ELSE LLRHO=.FALSE. END IF ELSE LLRHO=.FALSE. END IF C IF(LLRHO.OR.LLRLA.OR.LLRMU) THEN C Density and/or elastic perturbations are not zero C in this gridpoint. We will compute its contribution C to Born approximation. C C Loop over multivalued quantities of the Green function DO 29 JJ=IRAM(INUMG+I-1)+1,IRAM(INUMG+I) C Auxiliary names. TTG=RAM(IMTTG+JJ) C PG1=RAM(IPG1+JJ) PG2=RAM(IPG2+JJ) PG3=RAM(IPG3+JJ) C IF(L2D) THEN TTXXG=RAM(ITTXXG+JJ) END IF C GR11=RAM(IGR11+JJ) GR21=RAM(IGR21+JJ) GR31=RAM(IGR31+JJ) GR12=RAM(IGR12+JJ) GR22=RAM(IGR22+JJ) GR32=RAM(IGR32+JJ) GR13=RAM(IGR13+JJ) GR23=RAM(IGR23+JJ) GR33=RAM(IGR33+JJ) C GI11=RAM(IGI11+JJ) GI21=RAM(IGI21+JJ) GI31=RAM(IGI31+JJ) GI12=RAM(IGI12+JJ) GI22=RAM(IGI22+JJ) GI32=RAM(IGI32+JJ) GI13=RAM(IGI13+JJ) GI23=RAM(IGI23+JJ) GI33=RAM(IGI33+JJ) C Loop over multivalued quantities of the incident C wave field. DO 28 II=IRAM(INUMU+I-1)+1,IRAM(INUMU+I) C Auxiliary names. TTU=RAM(IMTTU+II) C PU1=RAM(IPU1+II) PU2=RAM(IPU2+II) PU3=RAM(IPU3+II) C IF(L2D) THEN TTXX=RAM(ITTXX+II) END IF C UR1=RAM(IUR1+II) UR2=RAM(IUR2+II) UR3=RAM(IUR3+II) C UI1=RAM(IUI1+II) UI2=RAM(IUI2+II) UI3=RAM(IUI3+II) C C Computing travel time from the source to the receiver. TT=TTU+TTG C And finding minimum and maximum travel time. TTMIN=AMIN1(TTMIN,TT) TTMAX=AMAX1(TTMAX,TT) C C Computation of auxiliary variables and computation of C frequency independent parts of the Born approximation: C B1R=0. B1I=0. B2R=0. B2I=0. B3R=0. B3I=0. C C ........................................................ C C 1-st term (and auxiliary variables for the 3-rd) IF(LLRHO.OR.LLRMU) THEN C Auxiliary variables for the 1-st and the 3-rd term AUAG1R=UR1*GR11+UR2*GR21+UR3*GR31 & -UI1*GI11-UI2*GI21-UI3*GI31 AUAG1I=UI1*GR11+UI2*GR21+UI3*GR31 & +UR1*GI11+UR2*GI21+UR3*GI31 C AUAG2R=UR1*GR12+UR2*GR22+UR3*GR32 & -UI1*GI12-UI2*GI22-UI3*GI32 AUAG2I=UI1*GR12+UI2*GR22+UI3*GR32 & +UR1*GI12+UR2*GI22+UR3*GI32 C AUAG3R=UR1*GR13+UR2*GR23+UR3*GR33 & -UI1*GI13-UI2*GI23-UI3*GI33 AUAG3I=UI1*GR13+UI2*GR23+UI3*GR33 & +UR1*GI13+UR2*GI23+UR3*GI33 C IF(LLRHO) THEN C Computation of the 1-st term B1R=RHO*AUAG1R B1I=RHO*AUAG1I B2R=RHO*AUAG2R B2I=RHO*AUAG2I B3R=RHO*AUAG3R B3I=RHO*AUAG3I END IF END IF C C ........................................................ C C 2-nd term IF(LLRLA) THEN C Auxiliary variables for the 2-nd term PGAG1R=PG1*GR11+PG2*GR21+PG3*GR31 PGAG1I=PG1*GI11+PG2*GI21+PG3*GI31 C PGAG2R=PG1*GR12+PG2*GR22+PG3*GR32 PGAG2I=PG1*GI12+PG2*GI22+PG3*GI32 C PGAG3R=PG1*GR13+PG2*GR23+PG3*GR33 PGAG3I=PG1*GI13+PG2*GI23+PG3*GI33 C PUAUR=PU1*UR1+PU2*UR2+PU3*UR3 PUAUI=PU1*UI1+PU2*UI2+PU3*UI3 C C Computation of the 2-st term B1R=B1R+RLA*(PGAG1R*PUAUR-PGAG1I*PUAUI) B1I=B1I+RLA*(PGAG1I*PUAUR+PGAG1R*PUAUI) B2R=B2R+RLA*(PGAG2R*PUAUR-PGAG2I*PUAUI) B2I=B2I+RLA*(PGAG2I*PUAUR+PGAG2R*PUAUI) B3R=B3R+RLA*(PGAG3R*PUAUR-PGAG3I*PUAUI) B3I=B3I+RLA*(PGAG3I*PUAUR+PGAG3R*PUAUI) END IF C C ........................................................ C C 3-rd term IF(LLRMU) THEN C Auxiliary variables for the 3-rd term PUAG1R=PU1*GR11+PU2*GR21+PU3*GR31 PUAG1I=PU1*GI11+PU2*GI21+PU3*GI31 C PUAG2R=PU1*GR12+PU2*GR22+PU3*GR32 PUAG2I=PU1*GI12+PU2*GI22+PU3*GI32 C PUAG3R=PU1*GR13+PU2*GR23+PU3*GR33 PUAG3I=PU1*GI13+PU2*GI23+PU3*GI33 C C PGAUR=PG1*UR1+PG2*UR2+PG3*UR3 PGAUI=PG1*UI1+PG2*UI2+PG3*UI3 C C PUPG=PU1*PG1+PU2*PG2+PU3*PG3 C C Computation of the 3-rd term B1R=B1R+RMU*(PUPG*AUAG1R & +PUAG1R*PGAUR-PUAG1I*PGAUI) B1I=B1I+RMU*(PUPG*AUAG1I & +PUAG1R*PGAUI+PUAG1I*PGAUR) B2R=B2R+RMU*(PUPG*AUAG2R & +PUAG2R*PGAUR-PUAG2I*PGAUI) B2I=B2I+RMU*(PUPG*AUAG2I & +PUAG2R*PGAUI+PUAG2I*PGAUR) B3R=B3R+RMU*(PUPG*AUAG3R & +PUAG3R*PGAUR-PUAG3I*PGAUI) B3I=B3I+RMU*(PUPG*AUAG3I & +PUAG3R*PGAUI+PUAG3I*PGAUR) END IF C C ........................................................ C C Computation of term C (frequency independent part)*exp(i*2*pi*fmin(tauU+tauG)) C Warning. DER here means real, DEI imaginary part of this C exponential. But DER and DEI will further have C a different meaning. DER=COS(P2I*FMIN*TT) DEI=SIN(P2I*FMIN*TT) IF(L2D) THEN DER=DER/SQRT(TTXX+TTXXG) DEI=DEI/SQRT(TTXX+TTXXG) END IF C BR=B1R*DER-B1I*DEI BI=B1R*DEI+B1I*DER B1R=BR B1I=BI C BR=B2R*DER-B2I*DEI BI=B2R*DEI+B2I*DER B2R=BR B2I=BI C BR=B3R*DER-B3I*DEI BI=B3R*DEI+B3I*DER B3R=BR B3I=BI C DER=COS(P2I*FSTEP*TT) DEI=SIN(P2I*FSTEP*TT) C The innermost loop over frequencies: DO 27 JF=IBORN+1,IBORN+N6F,6 C Computation of Born approximation (not multiplied C by term frequency squared and term D123) using C previously computed auxiliary variables: C Variables J2-J6 are auxiliary for better performance. J2=JF+1 J3=JF+2 J4=JF+3 J5=JF+4 J6=JF+5 C RAM(JF)=RAM(JF)+B1R RAM(J2)=RAM(J2)+B1I RAM(J3)=RAM(J3)+B2R RAM(J4)=RAM(J4)+B2I RAM(J5)=RAM(J5)+B3R RAM(J6)=RAM(J6)+B3I C BR=B1R*DER-B1I*DEI BI=B1R*DEI+B1I*DER B1R=BR B1I=BI C BR=B2R*DER-B2I*DEI BI=B2R*DEI+B2I*DER B2R=BR B2I=BI C BR=B3R*DER-B3I*DEI BI=B3R*DEI+B3I*DER B3R=BR B3I=BI 27 CONTINUE 28 CONTINUE 29 CONTINUE END IF 30 CONTINUE C C Adjustment of Born approximation. C Multiplication by term C 1D, 3D: Frequency squared and D1*D2*D3 (in the code denoted by C D123} C 2D: Frequency to the power of three halves, D1*D3 (in the code C denoted by D123), square root of pi and complex number 1+i. DO 31 JF=0,NF-1 F=FMIN+FLOAT(JF)*FSTEP C Variables J1-J6,J6F,DOM2 are auxiliary for better performance. IF(L2D) THEN DOM2=((SQRT(P2I*F))**3)*D123*SQRT(PI) ELSE DOM2=D123*(P2I*F)**2. END IF J6F=IBORN+6*JF J1 = 1 + J6F J2 = 2 + J6F J3 = 3 + J6F J4 = 4 + J6F J5 = 5 + J6F J6 = 6 + J6F C IF(L2D) THEN BR=DOM2*(RAM(J1)-RAM(J2)) BI=DOM2*(RAM(J1)+RAM(J2)) RAM(J1)=BR RAM(J2)=BI BR=DOM2*(RAM(J3)-RAM(J4)) BI=DOM2*(RAM(J3)+RAM(J4)) RAM(J3)=BR RAM(J4)=BI BR=DOM2*(RAM(J5)-RAM(J6)) BI=DOM2*(RAM(J5)+RAM(J6)) RAM(J5)=BR RAM(J6)=BI ELSE RAM(J1)=DOM2*RAM(J1) RAM(J2)=DOM2*RAM(J2) RAM(J3)=DOM2*RAM(J3) RAM(J4)=DOM2*RAM(J4) RAM(J5)=DOM2*RAM(J5) RAM(J6)=DOM2*RAM(J6) END IF 31 CONTINUE C C....................................................................... C C Saving results into the output file. C AMAX=0. DO 40 I=IBORN+1,IBORN+N6F,6 AMAX=AMAX1(ABS(RAM(I)) ,ABS(RAM(I+1)), * ABS(RAM(I+2)),ABS(RAM(I+3)), * ABS(RAM(I+4)),ABS(RAM(I+5)), * AMAX) 40 CONTINUE C IF(AMAX.LE.0.) THEN IF(TTMIN.GT.0.) THEN TTMAX=0. ELSE TTMAX=TTMIN-1. END IF END IF C IF(XR1.LT.-99.999.OR.999.999.LT.XR1.OR. * XR2.LT.-99.999.OR.999.999.LT.XR2.OR. * XR3.LT.-99.999.OR.999.999.LT.XR3) THEN WRITE(LU3,'(3A,3(F7.0,1X),3(E11.5,1X),A)') * '''',TXTREC,''' ',XR1,XR2,XR3,TTMIN,TTMAX,AMAX,'/' ELSE WRITE(LU3,'(3A,3(F7.3,1X),3(E11.5,1X),A)') * '''',TXTREC,''' ',XR1,XR2,XR3,TTMIN,TTMAX,AMAX,'/' END IF IF(TTMIN.LE.TTMAX) THEN DO 41 I=IBORN+1,IBORN+N6F,12 I01=IFIX(99999.1*RAM(I)/AMAX) I02=IFIX(99999.1*RAM(I+1)/AMAX) I03=IFIX(99999.1*RAM(I+2)/AMAX) I04=IFIX(99999.1*RAM(I+3)/AMAX) I05=IFIX(99999.1*RAM(I+4)/AMAX) I06=IFIX(99999.1*RAM(I+5)/AMAX) IF(I.LT.(IBORN+1+6*(NF-1))) THEN I07=IFIX(99999.1*RAM(I+6)/AMAX) I08=IFIX(99999.1*RAM(I+7)/AMAX) I09=IFIX(99999.1*RAM(I+8)/AMAX) I10=IFIX(99999.1*RAM(I+9)/AMAX) I11=IFIX(99999.1*RAM(I+10)/AMAX) I12=IFIX(99999.1*RAM(I+11)/AMAX) WRITE(LU3,'(12I6)') I01,I02,I03,I04,I05,I06, * I07,I08,I09,I10,I11,I12 ELSE WRITE(LU3,'(12I6)') I01,I02,I03,I04,I05,I06 END IF 41 CONTINUE END IF 50 CONTINUE C WRITE(LU3,'(A)') '/' C CLOSE(LU1) CLOSE(LU3) CLOSE(LU4) WRITE(*,'(A)') '+GRDBORN: Done. ' STOP END C C======================================================================= C SUBROUTINE CHLR(FNAME,TDEF,LU,N,LMSSG,ARRAY,INPUT) CHARACTER*80 MSSG CHARACTER*(*) FNAME,TDEF LOGICAL INPUT,LMSSG INTEGER LU,N,LENGTH REAL ARRAY(N) C C Subroutine CHeck and Load is designed to check whether there is any C input filename in the input SEP parameter file. Filename is specified C by parameter FNAME. If there is no input, variable INPUT will be false C otherwise true. C If variable LMSSG is true, INPUT.EQ..FALSE. creates error. If MSSG is C false, nothing will happen. C If there is any input filename, real values contained in this file C will be read from the file. Otherwise no reading will be performed. C C Input: C FNAME...String containing the name of the parameter. Except for C its case, it should match the parameter name in the input C SEP parameter file. C TDEF... Default value of the parameter. C LU... Logical unit number to be used. C N... Array dimension (number of elements to read). C LMSSG...Logical variable. Determines whether absence of input C file will create error or not. C C Output: C ARRAY.. Array having been read. C INPUT...Logical variable. TRUE if input file has been specified, C FALSE otherwise. C C----------------------------------------------------------------------- C C Other variables: CHARACTER*80 TOUT C TOUT... Value of the parameter specified by FNAME. C C----------------------------------------------------------------------- C INPUT=.TRUE. C CALL RSEP3T(FNAME,TOUT,TDEF) C IF(TOUT.EQ.' ') THEN INPUT=.FALSE. END IF C IF(INPUT) THEN CALL RARRAY(LU,TOUT,'FORMATTED',.TRUE.,0.,N,ARRAY) ELSE IF(LMSSG) THEN C GRDBORN-10 WRITE(MSSG,'(3A)') 'GRDBORN-10: Input file ', & FNAME(1:LENGTH(FNAME)), & ' not specified.' CALL ERROR(MSSG) END IF END IF END C C======================================================================= C SUBROUTINE CHLI(FNAME,TDEF,LU,N,LMSSG,IARRAY,INPUT) CHARACTER *80 MSSG CHARACTER*(*) FNAME,TDEF logical INPUT,LMSSG INTEGER LU,N,LENGTH INTEGER IARRAY(N) C C Subroutine CHeck and Load is designed to check whether there is any C input filename in the input SEP parameter file. Filename is specified C by parameter FNAME. If there is no input,variable INPUT will be false C otherwise true. C If variable LMSSG is true, INPUT.EQ..FALSE. creates error. If MSSG is C false, nothing will happen. C If there is any input filename, integer values contained in this file C will be read from the file. Otherwise no reading will be performed C C Input: C FNAME...String containing the name of the parameter. Except for C its case, it should match the parameter name in the input C SEP parameter file. C TDEF... Default value of the parameter. C LU... Logical unit number to be used. C N... Array dimension (number of elements to read). C LMSSG...Logical variable. Determines whether absence of input C file will create error or not. C C Output: C IARRAY..Array having been read. C INPUT...Logical variable. TRUE if input file has been specified, C FALSE otherwise. C C----------------------------------------------------------------------- C C Other variables: CHARACTER*80 TOUT C TOUT... Value of the parameter specified by FNAME. C C----------------------------------------------------------------------- C INPUT=.TRUE. C CALL RSEP3T(FNAME,TOUT,TDEF) C IF(TOUT.EQ.' ') THEN INPUT=.FALSE. END IF C IF(INPUT) THEN CALL RARRAI(LU,TOUT,'FORMATTED',.TRUE.,0,N,IARRAY) ELSE IF(LMSSG) THEN WRITE(MSSG,'(3A)') 'GRDBORN-10: Input file ', & FNAME(1:LENGTH(FNAME)), & ' not specified.' CALL ERROR(MSSG) END IF END IF END C C======================================================================= C SUBROUTINE CHLR2(FNAME,TOUT,LU,N,LMSSG,ARRAY,INPUT) CHARACTER *80 MSSG CHARACTER*(*) FNAME,TOUT LOGICAL INPUT,LMSSG INTEGER LU,N,LENGTH REAL ARRAY(N) C C Subroutine CHeck and Load is designed to check whether there is any C input filename. If not, variable INPUT will be false otherwise true. C If variable LMSSG is true, INPUT.EQ..FALSE. creates error. If MSSG is C false, nothing will happen. C If there is any input filename, real values contained in this file C will be read from the file. Otherwise no reading will be performed C C Input: C FNAME...String containing the name of the parameter specifying C the file C TOUT... String containg the filename. C LU... Logical unit number to be used. C N... Array dimension (number of elements to read). C LMSSG...Logical variable. Determines whether absence of input C file will create error or not. C C Output: C ARRAY.. Array having been read. C INPUT...Logical variable. TRUE if input file has been specified, C FALSE otherwise. C C----------------------------------------------------------------------- C INPUT=.TRUE. C IF(TOUT.EQ.' ') THEN INPUT=.FALSE. END IF C IF(INPUT) THEN CALL RARRAY(LU,TOUT,'FORMATTED',.TRUE.,0.,N,ARRAY) ELSE IF(LMSSG) THEN WRITE(MSSG,'(3A)') 'GRDBORN-10: Input file ', & FNAME(1:LENGTH(FNAME)), & ' not specified.' CALL ERROR(MSSG) END IF END IF END C C======================================================================= C SUBROUTINE CHLI2(FNAME,TOUT,LU,N,LMSSG,IARRAY,INPUT) CHARACTER *80 MSSG CHARACTER*(*) FNAME,TOUT LOGICAL INPUT,LMSSG INTEGER LU,N,LENGTH INTEGER IARRAY(N) C C Subroutine CHeck and Load is designed to check whether there is any C input filename. If not, variable INPUT will be false otherwise true. C If variable LMSSG is true, INPUT.EQ..FALSE. creates error. If MSSG is C false, nothing will happen. C If there is any input filename, integer values contained in this file C will be read from the file. Otherwise no reading will be performed C C Input: C FNAME...String containing the name of the parameter specifying C the file C TOUT... String containg the filename. C LU... Logical unit number to be used. C N... Array dimension (number of elements to read). C LMSSG...Logical variable. Determines whether absence of input file C will create error or not. C C Output: C IARRAY..Array having been read. C INPUT...Logical variable. TRUE if input file has been specified, C FALSE otherwise. C C----------------------------------------------------------------------- C INPUT=.TRUE. C IF(TOUT.EQ.' ') THEN INPUT=.FALSE. END IF C IF(INPUT) THEN CALL RARRAI(LU,TOUT,'FORMATTED',.TRUE.,0,N,IARRAY) ELSE IF(LMSSG) THEN WRITE(MSSG,'(3A)') 'GRDBORN-10: Input file ', & FNAME(1:LENGTH(FNAME)), & ' not specified.' CALL ERROR(MSSG) END IF END IF END C C======================================================================= C SUBROUTINE CHSP(MEMORY,NEEDED,QNAME) INTEGER MEMORY,NEEDED,LENGTH CHARACTER *(*) QNAME C C Subroutine designed to CHeck SPace in memory which we would like C to use to store data or as a computational space. C C Input: C MEMORY..Reserved memory space. C NEEDED..Needed memory space. C QNAME ..Name of the last quantity in the memory. C C----------------------------------------------------------------------- C C Other variables: CHARACTER *80 MMEM,MNEED CHARACTER *140 MSSG C C----------------------------------------------------------------------- C IF(MEMORY.LT.NEEDED) THEN WRITE(MMEM,*) MEMORY WRITE(MNEED,*) NEEDED C GRDBORN-11 WRITE(MSSG,'(7A)') 'GRDBORN-11: Not enough space for ', & QNAME(1:LENGTH(QNAME)),' in array RAM.', & ' Available:',MMEM(1:LENGTH(MMEM)), & ' Needed:',MNEED(1:LENGTH(MNEED)) CALL ERROR(MSSG) END IF END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for C C======================================================================= C