C
C Program INVPTS to calculate the derivatives of functions, describing C interfaces or material parameters, with respect to the model B-spline C coefficients C C Version: 5.50 C Date: 2001, June 9 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C E-mail: klimes@seis.karlov.mff.cuni.cz C C....................................................................... C C Program INVPTS assumes all model parameters (coefficients) stored in C the common block /VALC/ as in the submitted versions of user-defined C model specification FORTRAN77 source code files 'srfc.for', 'parm.for' C and 'val.for'. Thus, unlike the other parts of the complete ray C tracing, the INVTT program cannot work with user's modifications of C subroutines SRFC1, SRFC2, PARM1, and PARM2. C C The material parameters at the given points are converted into such C their powers which are interpolated by B-splines. Then the inversion C of the functions of coordinates is linear and one iteration yields C the exact solution within the rounding errors. On the other hand, C inversion of a material parameter depending on another material C parameters, e.g., RHO=W(VP(X1,X2,X3)), require the second iteration. C More than one iteration may also be required to reduce the rounding C errors. 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 Data specifying the model: C MODEL='string'... Name of the input data file describing the C model. For the description of the data refer to C subroutine file 'model.for'. C Default: MODEL='model.dat' C Data specifying other input files: C M2IN='string'... Name of the input file containing number M2IN of C values already processed by programs 'invpts.for' or C 'invtt.for', i.e., the name of output file M2 from the C previous execution of 'invpts.for' or 'invtt.for'. C The corresponding output values of this program will be C appended after M1*M2IN input values of file GM1 and after C M2IN input values of files GM2, GM3 and DM1, respectively. C Default: M2IN=' ' means M2IN=0 (no previous values). C PTS='string'... String with the name of the input data file C containing the individual given points with the values to C be fitted by the model. Points situated outside the model C box are not considered. C This program may be executed several times to take into C account several files with individual points, or the same C file with several values at each point (e.g., both P and C S wave velocities at each point). C If the filename is blank, no individual points are given. C Description of file PTS C Default: PTS=' ' C LIN='string'... Alternative to PTS, if the given points are C arranged into lines. Differs from PTS just by the C form of the file. C If the filename is blank, no lines of points are given. C Default: LIN=' ' C GRD='string'... String with the name of the input data file C containing the grided values to be fitted by the model. C The grid may also contain undefined values. C If the whole grid cell centred at a point is situated C outside the model box, the point is not considered. C This program may be executed several times to take into C account several grids of values. C If the filename is blank, no gridded values are assumed. C Default: GRD=' ' C GRDERR='string'... String with the name of the input data file C containing the grided standard deviations of the given C values. C If the filename is blank, unit standard deviations are C assumed. C The standard deviations are multiplied by the value of C ERRMUL, described later on. C Default: GRDERR=' ' C GRDICB='string'... String with the name of the input data file C containing the grided indices of the complex blocks C corresponding to file GRD, if file GRD contains the C material parameters of different complex blocks. C If GRDICB=' ', the index of the surface or complex block C corresponding to file GRD is common for all gridpoints C and is given by parameter INDFUN. C Default: GRDICB=' ' C For general description of the files with gridded data refer to C file forms.htm. C Data specifying output files: C M1='string'... Name of the output file containing number M1 of C model parameters (a single integer). The same file may be C generated by programs 'invsoft.for' and 'invtt.for'. C The file is not generated if the value of M1 is blank. C Default: M1=' ' C Note: Default of 'invsoft.for' is M1='m1.out', C default of 'invtt.for' is M1=' '. C M2='string'... Name of the output file containing number M2 of C accumulated values (a single integer). M2 is M2IN C increased by the number of given values. C The file is not generated if the value of M2 is blank. C Default: M2='m2.out' C Data specifying input/output files with matrix and vector components: C GM1='string'... String with the name of the input/output data file C to accumulate the matrix of the derivatives of M2 given C values with respect to M1 model coefficients. C The matrix has M1*M2IN components on input and M1*M2 C components on output. The formatted file is composed of C real-valued matrix components to be read at once by the C list-directed input. C If the filename is blank, no file is read nor written. C The file is not read, just written if M2IN=0. C Default: GM1=' ' C GM2='string'... String with the name of the input/output data file C containing the vector composed of differences between the C given values and the corresponding values calculated in C the model. The vector has M2IN components on input and M2 C components on output. The formatted file is composed C of real-valued vector components to be read at once by C list-directed input. C If the filename is blank, no file is read nor written. C The file is not read, just written if M2IN=0. C Default: GM2=' ' C GM3='string'... String with the name of the input/output data file C containing the vector composed of the values calculated in C the model. The vector has M2IN components on input and M2 C components on output. The formatted file is composed of C real-valued vector components to be read at once by C list-directed input. C If the filename is blank, no file is read nor written. C The file is not read, just written if M2IN=0. C Default: GM3=' ' C DM1='string'... String with the name of the input/output data file C containing the diagonal matrix composed of the variances C (squares of standard deviations) of the given values. C The diagonal has M2IN components on input and M2 C components on output. The formatted file is composed C of real-valued diagonal components to be read at once by C list-directed input. C If the filename is blank, no file is read nor written. C The file is not read, just written if M2IN=0. C Default: DM1=' ' C For general description of the files with matrices refer to file C forms.htm. C Form of the files with matrices: C FORMM='string' ... Form of the files with matrices. Allowed values C are FORMM='formatted' and FORMM='unformatted'. If the form C differs for input and for output files, FORMMR and FORMMW C should be used instead of FORMM. C Default: FORMM='formatted' C FORMMR='string' ... Form of the files with matrices to be read. C Default: FORMMR=FORMM C FORMMW='string' ... Form of the files with matrices to be written. C Default: FORMMW=FORMM C Data specifying the function corresponding to the given values: C ICLASS=integer... Class of model parameters to be inverted: C ICLASS=0: All model parameters are inverted. C ICLASS=1: Only model parameters describing interfaces are C inverted. C ICLASS=2: Only model parameters describing material C parameters are inverted. C Default: ICLASS=0 C MPAR=integer... Material parameter corresponding to the given C values. C MPAR=0: The given values correspond to the functions C describing the surfaces covering structural interfaces. C MPAR=positive integer: The given values describe a C material parameter: C MPAR=1: P wave velocity, C MPAR=2: S wave velocity, C MPAR=3: density, C MPAR=4: P wave loss factor, C MPAR=5: S wave loss factor, C MPAR=6 to 26: Reduced (i.e., divided by the density) C anisotropic elastic parameters A11, A12, A22, A13, C A23, A33, A14, A24, A34, A44, A15, A25, A35, A45, A55, C A16, A26, A36, A46, A56 or A66. C MPAR=27-47: Reduced (i.e., divided by the density) C imaginary parts of anisotropic elastic parameters Q11, C Q12, Q22, Q13, Q23, Q33, Q14, Q24, Q34, Q44, Q15, Q25, C Q35, Q45, Q55, Q16, Q26, Q36, Q46, Q56 or Q66. C Default: MPAR=0 C POWERM=real... The given values correspond to the POWERMth power C of the material parameter determined by MPAR. C Default: POWERM=1. C KOLFUN=integer... Specifies the column in input file 'PTS' or C 'LIN' containing the index of the surface (for MPAR=0) or C the complex block. C If KOLFUN=0, the index is common for all points and is C given by parameter INDFUN. C Default: KOLFUN=0 C INDFUN=integer... Index of the surface (for MPAR=0) or of the C complex block to which the given values correspond. C Must be specified and positive if KOLFUN=0 for individual C points or GRDICB=' ' for gridded values. C Default: INDFUN=0 C Data specifying the values to be processed and their accuracy: C KOLUMN=integer... Specifies the column in input file 'PTS' or C 'LIN' containing the given values. Note that the first C 3 columns contain coordinates of the points. C If KOLUMN=0, zero values are assumed (option often used C for the points at interfaces, where the function C describing the corresponding surface should be zero). C Default: KOLUMN=0 C KOLERR=integer... Specifies the column in input file 'PTS' or C 'LIN' containing the standard deviations of the given C values. C If KOLERR=0, unit standard deviations are assumed. C Note that the standard deviations are multiplied by C ERRMUL, described below. C Default: KOLERR=0 C ERRMUL=real... Multiplication factor for the standard deviations C of the given values. Often used with GRDERR=' ' to C specify the standard deviation constant over the grid or C with KOLERR=0 to specify the standard deviation common for C all points. C Default: ERRMUL=1. C C C Input file PTS with the points: C (1) None to several strings terminated by / (a slash) C (2) For each point data: C (2.1) 'NAME',X1,X2,X3,V1,...,VN,/ C 'NAME'... Name of the point. Not considered. May be blank. C X1,X2,X3... Coordinates of the point C V1,...,VN...Optional given values and their standard deviations, C see input parameters KOLUMN, KOLERR and KOLFUN. C /... Values must be terminated by a slash. C (3) / or end of file. C C C Input file LIN with the lines: C (1) None to several strings terminated by / (a slash) C (2) For each line data (2.1), (2.2) and (2.3): C (2.1) 'NAME',X1,X2,X3,/ C 'NAME'... Name of the line. Not considered. May be blank but C must be different from '$'. C X1,X2,X3... Optional coordinates of the reference point of the C line. Not considered. C /... List of values must be terminated by a slash. C (2.2) For each point of the line data (2.2.1): C (2.2.1) X1,X2,X3,V1,...,VN,/ C X1,X2,X3... Coordinates of the point of the line. C V1,...,VN...Optional given values and their standard deviations, C see input parameters KOLUMN, KOLERR and KOLFUN. C /... List of values must be terminated by a slash. C (2.3) / C (3) / or end of file. C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C C Common blocks /MODELC/ and /VALC/: INCLUDE 'model.inc' C model.inc INCLUDE 'val.inc' C val.inc C None of the storage locations of the common block are altered. C C External procedures directly referred: EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,RARRAI,RARRAY,WARRAY EXTERNAL LENGTH,MODEL1,SRFC2,PARM2,PARM3,SOFT,VAR6,RMAT,WMAT INTEGER LENGTH C C....................................................................... C C Constants: INTEGER LU1 PARAMETER (LU1=11) REAL UNDEF PARAMETER (UNDEF=-9.999E9) C C Filenames: CHARACTER*80 FILE1 C C Data: INTEGER N1,N2,N3,M1,M2,M2IN INTEGER ICLASS,MPAR,KOLUMN,KOLERR,KOLFUN,INDFUN REAL D1,D2,D3,O1,O2,O3,POWERM,ERRMUL C C Other variables: CHARACTER*1 TEXT INTEGER N1N2N3,KOLMAX,NFREE,NUNDEF,NOFF,NPTS6,NFUN INTEGER I,I1,I2,I3,J1,J2,J3,K1,K2,K3 REAL F(10,47),W(47,2),POWER(47),X1,X2,X3,B0I,AUX,CS(1) EQUIVALENCE (F(1,1),W(1,1)) C C Usage of RAM: C RAM(1:6*NPTS): C RAM(1,*)=X1 C RAM(2,*)=X2 C RAM(3,*)=X3 C RAM(4,*)=given value, later given minus reference value C IRAM(4,*)=MPAR, later RAM(4,*)=reference value C RAM(6,*)=standard deviation, later the variance C RAM(6*NPTS+1:6*NPTS+M1*M2): Derivatives C RAM(6*NPTS+M1*M2+1:6*NPTS+M1*M2+M1): Indices of model coefficients C C----------------------------------------------------------------------- C C Opening data files and reading the input data: C C Reading main input data: WRITE(*,'(A)') '+INVPTS: Enter input filename: ' FILE1=' ' READ (*,*) FILE1 IF(FILE1.EQ.' ') THEN C INVPTS-01 CALL ERROR('INVPTS-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. END IF WRITE(*,'(A)') '+INVPTS: Working... ' CALL RSEP1(LU1,FILE1) C C Reading input data for the model: CALL RSEP3T('MODEL' ,FILE1 ,'model.dat') OPEN(LU1,FILE=FILE1,STATUS='OLD') CALL MODEL1(LU1) CLOSE(LU1) C C Reading input parameters describing the input values: CALL RSEP3I('KOLUMN',KOLUMN,0) CALL RSEP3I('KOLERR',KOLERR,0) CALL RSEP3I('KOLFUN',KOLFUN,0) CALL RSEP3I('INDFUN',INDFUN,0) KOLMAX=MAX0(3,KOLUMN,KOLERR,KOLFUN) C C....................................................................... C C Reading the points: C C Number of valid points NPTS6=0 C Number of points with zero indeces (e.g., situated in free space) NFREE=0 C Number of points with positive indices but undefined values NUNDEF=0 C Number of defined points with positive indices outside the model NOFF=0 C C Reading gridded values: CALL RSEP3T('GRD',FILE1 ,' ') IF(FILE1.NE.' ') THEN CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3R('D1',D1,1.) CALL RSEP3R('D2',D2,1.) CALL RSEP3R('D3',D3,1.) CALL RSEP3R('O1',O1,0.) CALL RSEP3R('O2',O2,0.) CALL RSEP3R('O3',O3,0.) N1N2N3=N1*N2*N3 IF(6*N1N2N3.GT.MRAM) THEN GO TO 99 END IF CALL RARRAY(LU1,FILE1,'FORMATTED',.TRUE.,UNDEF,N1N2N3,RAM(1)) C Reading indices of surfaces or complex blocks CALL RSEP3T('GRDICB',FILE1 ,' ') IF(FILE1.EQ.' '.AND.INDFUN.LE.0) THEN C INVPTS-03 CALL ERROR('INVPTS-03: Neither GRDICB nor INDFUN specified') C See the description of the input data. END IF IF(FILE1.EQ.' ') THEN DO 11 I=N1N2N3+1,2*N1N2N3 IRAM(I)=INDFUN 11 CONTINUE ELSE CALL RARRAI(LU1,FILE1,'FORMATTED',.TRUE.,0, * N1N2N3,IRAM(N1N2N3+1)) END IF C Reading errors CALL RSEP3T('GRDERR',FILE1 ,' ') IF(FILE1.EQ.' ') THEN DO 12 I=2*N1N2N3+1,3*N1N2N3 RAM(I)=1. 12 CONTINUE ELSE CALL RARRAY(LU1,FILE1,'FORMATTED',.TRUE.,1., * N1N2N3,RAM(2*N1N2N3+1)) END IF C Rearranging the values in the memory DO 13 I=N1N2N3+1,2*N1N2N3 RAM (3*I-2)=RAM (I-N1N2N3) IRAM(3*I-1)=IRAM(I ) RAM (3*I )=RAM (I+N1N2N3) 13 CONTINUE C Removing undefined values and values outside the model J1=NINT((BOUNDM(1)-O1)/D1) K1=NINT((BOUNDM(2)-O1)/D1) J2=NINT((BOUNDM(3)-O2)/D2) K2=NINT((BOUNDM(4)-O2)/D2) J3=NINT((BOUNDM(5)-O3)/D3) K3=NINT((BOUNDM(6)-O3)/D3) I=3*N1N2N3 DO 23 I3=0,N3-1 X3=O3+D3*FLOAT(I3) DO 22 I2=0,N2-1 X2=O2+D2*FLOAT(I2) DO 21 I1=0,N1-1 X1=O1+D1*FLOAT(I1) I=I+3 IF(IRAM(I-1).LE.0) THEN NFREE=NFREE+1 ELSE IF(RAM(I-2).EQ.UNDEF) THEN NUNDEF=NUNDEF+1 ELSE IF(J3.LE.I3.AND.I3.LE.K3) THEN IF(J2.LE.I2.AND.I2.LE.K2) THEN IF(J1.LE.I1.AND.I1.LE.K1) THEN RAM (NPTS6+1)=X1 RAM (NPTS6+2)=X2 RAM (NPTS6+3)=X3 RAM (NPTS6+4)=RAM (I-2) IRAM(NPTS6+5)=IRAM(I-1) RAM (NPTS6+6)=RAM (I ) NPTS6=NPTS6+6 END IF END IF END IF END IF 21 CONTINUE 22 CONTINUE 23 CONTINUE NOFF=N1N2N3-NFREE-NUNDEF-NPTS6/6 END IF C C Reading individual points: CALL RSEP3T('PTS',FILE1 ,' ') IF(FILE1.NE.' ') THEN IF(KOLFUN.LE.0.AND.INDFUN.LE.0) THEN C INVPTS-04 CALL ERROR('INVPTS-04: Neither KOLFUN nor INDFUN specified') C See the description of the input data. END IF OPEN(LU1,FILE=FILE1,STATUS='OLD') READ(LU1,*,END=39) (TEXT,I=1,20) C Loop over points 31 CONTINUE IF(NPTS6+3+KOLMAX.GT.MRAM) THEN GO TO 99 END IF TEXT='$' X1=0. X2=0. X3=0. READ(LU1,*,END=39) TEXT,X1,X2,X3, * (RAM(I),I=NPTS6+7,NPTS6+3+KOLMAX) IF(TEXT.EQ.'$') THEN GO TO 39 END IF NOFF=NOFF+1 IF(BOUNDM(1).LE.X1.AND.X1.LE.BOUNDM(2)) THEN IF(BOUNDM(3).LE.X2.AND.X2.LE.BOUNDM(4)) THEN IF(BOUNDM(5).LE.X3.AND.X3.LE.BOUNDM(6)) THEN RAM(NPTS6+1)=X1 RAM(NPTS6+2)=X2 RAM(NPTS6+3)=X3 IF(KOLUMN.GT.0) THEN RAM(NPTS6+4)=RAM(NPTS6+3+KOLUMN) ELSE RAM(NPTS6+4)=0. END IF IF(KOLFUN.GT.0) THEN IRAM(NPTS6+5)=NINT(RAM(NPTS6+3+KOLFUN)) ELSE IRAM(NPTS6+5)=INDFUN END IF IF(KOLERR.GT.0) THEN RAM(NPTS6+6)=RAM(NPTS6+3+KOLERR) ELSE RAM(NPTS6+6)=1. END IF NPTS6=NPTS6+6 NOFF=NOFF-1 END IF END IF END IF GO TO 31 C End of the loop over points 39 CONTINUE CLOSE(LU1) END IF C C Reading the lines of points: CALL RSEP3T('LIN',FILE1 ,' ') IF(FILE1.NE.' ') THEN IF(KOLFUN.LE.0.AND.INDFUN.LE.0) THEN C INVPTS-05 CALL ERROR('INVPTS-05: Neither KOLFUN nor INDFUN specified') C See the description of the input data. END IF OPEN(LU1,FILE=FILE1,STATUS='OLD') READ(LU1,*,END=49) (TEXT,I=1,20) C Loop over lines 41 CONTINUE TEXT='$' READ(LU1,*,END=90) TEXT,AUX,AUX,AUX IF(TEXT.EQ.'$') THEN GO TO 49 END IF C Loop over points 42 CONTINUE IF(NPTS6+3+KOLMAX.GT.MRAM) THEN GO TO 99 END IF X1=UNDEF X2=0. X3=0. READ(LU1,*,END=49) TEXT,X1,X2,X3, * (RAM(I),I=NPTS6+7,NPTS6+3+KOLMAX) IF(X1.EQ.UNDEF) THEN C End of the line GO TO 48 END IF NOFF=NOFF+1 IF(BOUNDM(1).LE.X1.AND.X1.LE.BOUNDM(2)) THEN IF(BOUNDM(3).LE.X2.AND.X2.LE.BOUNDM(4)) THEN IF(BOUNDM(5).LE.X3.AND.X3.LE.BOUNDM(6)) THEN RAM(NPTS6+1)=X1 RAM(NPTS6+2)=X2 RAM(NPTS6+3)=X3 IF(KOLUMN.GT.0) THEN RAM(NPTS6+4)=RAM(NPTS6+3+KOLUMN) ELSE RAM(NPTS6+4)=0. END IF IF(KOLFUN.GT.0) THEN IRAM(NPTS6+5)=NINT(RAM(NPTS6+3+KOLFUN)) ELSE IRAM(NPTS6+5)=INDFUN END IF IF(KOLERR.GT.0) THEN RAM(NPTS6+6)=RAM(NPTS6+3+KOLERR) ELSE RAM(NPTS6+6)=1. END IF NPTS6=NPTS6+6 NOFF=NOFF-1 END IF END IF END IF GO TO 42 C End of the loop over points 48 CONTINUE GO TO 41 C End of the loop over lines 49 CONTINUE CLOSE(LU1) END IF C C....................................................................... C C Matrix dimensions: C C Reading number of points processed previously: M2IN=0 CALL RSEP3T('M2IN',FILE1,' ') IF(FILE1.NE.' ') THEN OPEN(LU1,FILE=FILE1,STATUS='OLD') READ(LU1,*) M2IN CLOSE(LU1) END IF C C Writing the total number of points: M2=M2IN+NPTS6/6 CALL RSEP3T('M2',FILE1,'m2.out') IF(FILE1.NE.' ') THEN OPEN(LU1,FILE=FILE1) WRITE(LU1,'(I10)') M2 CLOSE(LU1) END IF C C Number of unknown model parameters: CALL RSEP3I('ICLASS',ICLASS,0) IF(ICLASS.LT.0.OR.2.LT.ICLASS) THEN C INVPTS-06 CALL ERROR('INVPTS-06: Incorrect class index ICLASS') C The value of ICLASS must be 0, 1 or 2. C Check the input data. END IF DO 51 I=1,47 W(I,1)=0. W(I,2)=0. 51 CONTINUE M1=0 IF(ICLASS.LE.1) THEN CALL SOFT(1,0,0,0,0,0,0,47,W,M1,IRAM(NPTS6+1),CS(1),1,CS(1)) END IF IF(ICLASS.EQ.0.OR.ICLASS.EQ.2) THEN CALL SOFT(2,0,0,0,0,0,0,47,W,M1,IRAM(NPTS6+1),CS(1),1,CS(1)) END IF IF(NPTS6+M1*M2+M1.GT.MRAM) THEN GO TO 99 END IF I1=0 IF(ICLASS.LE.1) THEN CALL SOFT(1,0,0,0,0,0,0,47,W,I1,IRAM(NPTS6+M1*M2+1), * CS(1),1,CS(1)) END IF IF(ICLASS.EQ.0.OR.ICLASS.EQ.2) THEN CALL SOFT(2,0,0,0,0,0,0,47,W,I1,IRAM(NPTS6+M1*M2+1), * CS(1),1,CS(1)) END IF CALL RSEP3T('M1',FILE1,' ') IF(FILE1.NE.' ') THEN OPEN(LU1,FILE=FILE1) WRITE(LU1,'(I10)') M1 CLOSE(LU1) END IF C C Information written to the screen: IF(NFREE.GT.0) THEN WRITE(*,'(A,I8,A)')'+INVPTS:',NFREE,' gridpoints in free space,' WRITE(*,'(A)') ' INVPTS: Working...' END IF IF(NUNDEF.GT.0) THEN WRITE(*,'(A,I8,A)') '+INVPTS:',NUNDEF,' undefined grid values,' WRITE(*,'(A)') ' INVPTS: Working...' END IF IF(NOFF.GT.0) THEN WRITE(*,'(A,I8,A)') '+INVPTS:',NOFF,' points outside model,' WRITE(*,'(A)') ' INVPTS: Working...' END IF C C....................................................................... C C Calculating the functional values and the derivatives with respect C to the model coefficients at the stored points: C CALL RSEP3I('MPAR' ,MPAR ,0) IF(MPAR.LT.0.OR.47.LT.MPAR) THEN C INVPTS-07 CALL ERROR('INVPTS-07: Wrong material parameter') C Material parameters are indexed by integers from 1 to 47, C but index greater than 47 has been encountered. C Check the input data. END IF CALL RSEP3R('POWERM',POWERM,1.) CALL RSEP3R('ERRMUL',ERRMUL,1.) C C Loop over the points: DO 69 I3=0,NPTS6-6,6 IF((ICLASS.LE.1.AND.MPAR.LE.0).OR. * ((ICLASS.EQ.0.OR.ICLASS.EQ.2).AND.MPAR.GE.1)) THEN C Evaluation of the function: IF(MPAR.LE.0) THEN C Functions describing surfaces: C CALL SRFC2(IRAM(I3+5),RAM(I3+1),F) CALL VAL2(1,IRAM(I3+5),1,RAM(I3+1),F,POWER) RAM(I3+4)=RAM(I3+4)-F(1,1) RAM(I3+5)=F(1,1) ELSE C Material parameters: CALL VAL2(2,IRAM(I3+5),47,RAM(I3+1),F,POWER) AUX=POWER(MPAR)/POWERM IF(AUX.NE.1.) THEN C Power of the given value RAM(I3+4)=RAM(I3+4)**AUX C Standard deviation of the power RAM(I3+6)=RAM(I3+6)*AUX*RAM(I3+4)**(AUX-1.) END IF RAM(I3+4)=RAM(I3+4)-F(1,MPAR) RAM(I3+5)=F(1,MPAR) END IF C C Variance (square of the standard deviation): RAM(I3+6)=(RAM(I3+6)*ERRMUL)**2 C C Derivatives of the function with respect to model coefficients: DO 61 I1=NPTS6+M1*(M2IN+I3/6)+1,NPTS6+M1*(M2IN+I3/6)+M1 RAM(I1)=0. 61 CONTINUE I2=0 62 CONTINUE I2=I2+1 CALL VAR6(MAX0(1,MPAR),I2,NFUN,I,B0I,AUX,AUX,AUX) IF(I2.LE.NFUN) THEN DO 63 I1=NPTS6+M1*M2+1,NPTS6+M1*M2+M1 IF(I.EQ.IRAM(I1)) THEN RAM(M1*(I3-NPTS6)/6+I1)=B0I GO TO 64 END IF 63 CONTINUE C INVPTS-08 CALL ERROR * ('INVPTS-08: Incorrect index of model parameter') 64 CONTINUE END IF IF(I2.LT.NFUN) GO TO 62 C END IF 69 CONTINUE C C....................................................................... C C Writing output files: C C Writing the derivatives: CALL RSEP3T('GM1',FILE1 ,' ') IF(FILE1.NE.' ') THEN IF (M2IN.GT.0) THEN CALL RMAT(LU1,FILE1,M1,M2IN,RAM(NPTS6+1)) END IF CALL WMAT(LU1,FILE1,M1,M2,RAM(NPTS6+1)) END IF C C Writing the given values minus the reference values: CALL RSEP3T('GM2',FILE1 ,' ') IF(FILE1.NE.' ') THEN IF(M2IN.GT.0) THEN CALL RMAT(LU1,FILE1,M2IN,1,RAM(NPTS6+1)) END IF DO 81 I=1,NPTS6/6 RAM(NPTS6+M2IN+I)=RAM(6*I-2) 81 CONTINUE CALL WMAT(LU1,FILE1,M2,1,RAM(NPTS6+1)) END IF C C Writing the reference values: CALL RSEP3T('GM3',FILE1 ,' ') IF(FILE1.NE.' ') THEN IF(M2IN.GT.0) THEN CALL RMAT(LU1,FILE1,M2IN,1,RAM(NPTS6+1)) END IF DO 82 I=1,NPTS6/6 RAM(NPTS6+M2IN+I)=RAM(6*I-1) 82 CONTINUE CALL WMAT(LU1,FILE1,M2,1,RAM(NPTS6+1)) END IF C C Writing the variances: CALL RSEP3T('DM1',FILE1 ,' ') IF(FILE1.NE.' ') THEN IF(M2IN.GT.0) THEN CALL RMAT(LU1,FILE1,M2IN,1,RAM(NPTS6+1)) END IF DO 83 I=1,NPTS6/6 RAM(NPTS6+M2IN+I)=RAM(6*I) 83 CONTINUE CALL WMAT(LU1,FILE1,M2,1,RAM(NPTS6+1)) END IF C 90 CONTINUE WRITE(*,'(A,I8,A)') '+INVPTS:',NPTS6/6,' points processed.' STOP C 99 CONTINUE C INVPTS-02 CALL ERROR('INVPTS-02: Too small dimension MRAM of array RAM') C Dimension MRAM of array RAM in include file C ram.inc should be increased. C MRAM should be MAX(6*N1*N2*N3,6*NPTS+M1*M2+M1) at the least. C Here NPTS=M2-M2IN is the number of valid points, i.e., of points C situated inside the model box, with defined given values. 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 INCLUDE 'modelv.for' C modelv.for INCLUDE 'metric.for' C metric.for INCLUDE 'srfc.for' C srfc.for INCLUDE 'parmv.for' C parmv.for INCLUDE 'valv.for' C valv.for INCLUDE 'fitv.for' C fitv.for INCLUDE 'var.for' C var.for INCLUDE 'soft.for' C soft.for INCLUDE 'spsp.for' C spsp.for C C======================================================================= C