C Subroutine file 'val.for' for function specification and interpolation C - designed to perform the interpolation of a set of functions in a C rectangular grid. C C Date: 1999, September 15 C Coded by Ludek Klimes C C....................................................................... C C This file consists of the following subroutines: C VALB... Block data subroutine initiating common block /VALC/ to C store the data describing the interpolated functions. C VALB C VAL1... Subroutine designed to read the input data for the C functions, to compute the coefficients of the expansion C and to store them in the memory. C VAL1 C SORTV...Auxiliary subroutine to subroutine VAL1. C SORTV C READV...Auxiliary subroutine to subroutine VAL1. C READV C VAL2... Subroutine evaluating the functions including their first C and second derivatives. The functions may be used to C specify various surfaces in the model, the space C distributions of various parameters, e.t.c. The functions C are represented as a tensor product of splines under C tension of at most 3 independent variables (i.e. a linear C combination of products of B-splines under tension of at C most 3 independent variables). Outside the specified C rectangular grid, the functions are extrapolated by their C analytic continuation. See lines of subroutine VAL2 with C 'CV3' in first 3 columns for comparison with the kind of C extrapolation used in old versions (outside the specified C rectangular grid, the functions were linear along straight C lines perpendicular to the boundary of the grid). The C functions may be embedded: the independent variable of C the function may be another function of the same group C (see below) foregoing in the input data. C VAL2 C Subroutines of this file employ routines CURVN1 (or its alternative C CURVB1), CURV2D (or its alternative CURVBD), SURFB1, SURFBD, VAL3B1, C VAL3BD, VGEN, TERMS, SNHCSH, TRIDEC, TRISOL, DSPLNZ, INTRVL from the C subroutine package 'FITPACK' by Alan Kaylor Cline, Department of C Computer Sciences, University of Texas at Austin. C C Note: C The lines denoted by '*V' in the first two columns of file C 'val.for' in subroutine VAL2 are designed to calculate the model C variations with respect to the model parameters. C File 'valv.for', intended for the model inversion, is created C from 'val.for' by replacing each '*V' in the first two columns C by spaces using program 'clean.for'. Subroutines VAR4 and VAR5 C of file 'var.for' may then be called to handle the variations. C C....................................................................... C C Classes of functions: C The interpolated functions are divided into some classes, e.g., C functions describing interfaces, functions describing the medium C parameters, functions describing the properties of the source, C etc. The number MCLASS of the defined classes is stored in the C memory and is initially zero. The new class may be defined by C means of the invocation of subroutine VAL1, see its input argument C ICLASS. During one invocation of VAL1, only the groups of C functions relevant to one class may be defined. Subroutine VAL1 C may be called several times even for one class to define its C groups successively, stage by stage. In this case, the input data C for the groups of functions relevant to one class may be read in C from various files. C C Groups of functions: C The interpolated functions of each class are divided into some C groups. For instance, the class of functions describing the medium C parameters is divided into groups corresponding to individual C complex blocks. The individual groups need not contain the same C number of functions. The group corresponding to a complex block C may contain e.g. the functions describing P-wave velocity, S-wave C velocity, density, etc. The functions not specified by the input C data but required by the program are defined and are zero. C C C Input data (read in by subroutine VAL1): C These input data define the groups of functions from the specified C class. The index ICLASS of the class is an input argument of C subroutine VAL1. If the class is not defined by a previous C invocation of VAL1, it is created. The number NGROUP of the C groups to be defined is an input argument of subroutine VAL1. The C data are read in by the list directed input (free format). C (1) NGROUP-times (i.e. once for each group of functions) input data C (1A)+(1B): C (1A) TEXTG,IGROUP C Identification of the group. C TEXTG...Any string. Its first 3 characters must differ from 'END' C and from any string identifying a physical quantity, C defined by input array TFUNCT of subroutine VAL1 (see C below). C IGROUP..Sequential number of the group in the class. C (1B) Several times Input data for one function, C see below. C If input array TFUNCT of subroutine VAL1 is fully filled by C spaces, 'Input data for one function' must be included C just NFUNCT-times (NFUNCT is an input argument of C subroutine VAL1). In this case, the input functions are C not identified by a string (see (1) of 'Input data for one C function'), their number and order must be a priori known. C This is, for instance, the case of smooth surfaces: each C group corresponds to one surface and contains just one C function. The index of the surface coincides with the C index of the group, and no identification and sorting of C functions inside a group is needed. C If input array TFUNCT of subroutine VAL1 is not fully filled by C spaces, 'Input data for one function' may be included C N-times, where 0.LE.N.LE.NFUNCT. In this case, the input C functions of individual groups are identified by a string C (see (1) of 'Input data for one function') in the input C data, their number and order may be arbitrary (note that C their number must be less than or equal to NFUNCT). This C is, for instance, the case of complex blocks: each group C corresponds to one complex block and contains some number C of functions describing material parameters. Individual C functions (material parameters) are identified by a string C in the input data. C (2) TEXTE,AUX C End of data. C TEXTE...String, the first 3 characters of which must be upper-case C 'END'. C AUX... Any number or a slash. C C C Input data for one function: C The data are read in by the list directed input (free format). In C the list of input data below, each numbered paragraph indicates C the beginning of a new input operation (new READ statement). If C the first letter of the symbolic name of the input variable is C I-N, the corresponding value in input data must be of the type C INTEGER. Otherwise (except TEXTF), the input parameter is of the C type REAL. C (1) TEXTF,POWER C Physical meaning of the function. This input is not performed if C input character array TFUNCT of subroutine VAL1 (see below) is C fully filled by spaces. C TEXTF...String identifying which physical quantity the function C describes. Only the first 3 characters are significant. C The first 3 characters of the string must not be equal to C 'END'. The set of meaningful strings is defined by input C array TFUNCT of subroutine VAL1 (see below). C POWER...The specified function is equal to the POWER-th power of C the physical quantity. C Default: POWER=1. C (2) IVAR1,IVAR2,IVAR3,SIGMA,POWERW,/ C The form of the function. C IVAR1,IVAR2,IVAR3... Denote the form of the function. The function C must be of the form C F(X1,X2,X3) =W(A1,A2,A3)-B1-B2-B3 . C X1, X2, X3 are the general coordinates. Each of A1, A2, C A3, B1, B2, B3 must be either: (a) one of general C coordinates X1, X2, X3, (b) another previously defined C function F(X1,X2,X3) of the same group, or (c) must be C left out. At most 3 of parameters A1-B3 may be of kind C (a) or (b). Note that IVAR1 controls the type of A1 and C B1, IVAR2 controls the type of A2 and B2, IVAR3 controls C the type of A3 and B3. C For IVAR1.EQ.0: A1, B1 are empty (left out), C for IVAR1.EQ.1: A1=X1, B1 is empty, C for IVAR1.EQ.2: A1=X2, B1 is empty, C for IVAR1.EQ.3: A1=X3, B1 is empty, C for IVAR1.GE.4: A1=F(X1,X2,X3), where F(X1,X2,X3) is C another function of the same group defined in the input C data as the (IVAR1-3)-th function of the group. B1 is C empty, C for IVAR1.EQ.-1: B1=X1, A1 is empty, C for IVAR1.EQ.-2: B1=X2, A1 is empty, C for IVAR1.EQ.-3: B1=X3, A1 is empty, C for IVAR1.LE.-4: B1=F(X1,X2,X3), where F(X1,X2,X3) is C another function of the same group defined in the input C data as the (-IVAR1-3)-th function of the group. A1 is C empty. C The meaning of the parameters IVAR2, IVAR3 is similar. C Examples: C IVAR1: IVAR2: IVAR3: the form of the function: C 1 2 3 F(X1,X2,X3)=W(X1,X2,X3) C 3 1 2 F(X1,X2,X3)=W(X3,X1,X2) C 1 2 0 F(X1,X2,X3)=W(X1,X2) C 1 2 -3 F(X1,X2,X3)=W(X1,X2)-X3 C 1 -3 2 F(X1,X2,X3)=W(X1,X2)-X3 C 5 0 0 F(X1,X2,X3)=W(F2(X1,X2,X3)), where C F2(X1,X2,X3) is the second function of the group defined C in the input data. Function W is interpolated by means of C splines under tension. C SIGMA...Is the tension factor (its sign is ignored). This value C indicates the curviness desired. If ABS(SIGMA) is nearly C zero (e.g. 0.001), the resulting surface is approximately C the tensor product of cubic splines. If ABS(SIGMA) is C large (e.g. 50.), the resulting surface is approximately C tri-linear. If SIGMA equals zero, tensor products of C cubic splines result. A recommended value for SIGMA is C approximately 1. In absolute value. C POWERW..Given grid values (7) correspond to the POWERW-th power of C interpolated function W. The given grid values (7) are C thus raised to the (1/POWERW)-th power immediately after C reading and then interpolated. C /... Obligatory slash at the end of line for future extensions. C Default: IVAR1=0, IVAR2=0, IVAR3=0, SIGMA=0, POWERW=1. C (3) NX(1),...,NX(NVAR) C The numbers of grid coordinates for the interpolation. C This input is performed if at least one of IVAR1, IVAR2, IVAR3 is C positive. C Each of NX(1),...,NX(NVAR) corresponds to one positive value of C IVAR1, IVAR2, IVAR3 and specifies the number of grid coordinates C corresponding to that independent variable of function W, see (2). C The sign of NX(1),...,NX(NVAR) is ignored. NVAR (.LE.3) is the C number of positive values of the above quantities IVAR1, IVAR2, C IVAR3, i.e. the number of independent variables of function W, C see (2). C (4) X1(1),...,X1(NX(1)) C The grid coordinates corresponding to the first independent C variable of function W, see (2). C This input is performed if NX(1) is specified, see (3), and is not C zero. The grid coordinates may be specified in any order. C (5) X2(1),...,X2(NX(2)) C The grid coordinates corresponding to the second independent C variable of function W, see (2). C This input is performed if NX(2) is specified, see (3), and is not C zero. The grid coordinates may be specified in any order. C (6) X3(1),...,X3(NX(3)) C The grid coordinates corresponding to the third independent C variable of function W, see (2). C This input is performed if NX(3) is specified, see (3), and is not C zero. The grid coordinates may be specified in any order. C (7) (((W(I,J,K),I=1,MAX(NX(1),1)),J=1,MAX(NX(2),1)),K=1,MAX(NX(3),1)) C The values of the POWERW-th power of function W at grid points. C Function value W(I,J,K) corresponds to point (X1(I),X2(J),X3(K)). C C======================================================================= C C C C Storage in the memory: C The parameters describing the interpolated functions are stored C in common block /VALC/ initialized in the following subroutine: C ------------------------------------------------------------------ BLOCK DATA VALB INCLUDE 'val.inc' C val.inc DATA IPAR(0)/0/ END C ------------------------------------------------------------------ C C======================================================================= C C C SUBROUTINE VAL1(LUN,ICLASS,NGROUP,NFUNCT,TFUNCT) INTEGER LUN,ICLASS,NGROUP,NFUNCT CHARACTER*(*) TFUNCT(NFUNCT) C C This subroutine reads the input data for a set of functions, C determines the parameters necessary to compute an interpolatory C function on a three-dimensional rectangular grid, and stores them in C the memory. The function determined can be represented as a tensor C product of splines under tension of at most 3 independent variables C (i.e. a linear combination of products of B-splines under tension of C at most 3 independent variables). The functions may be embedded: the C independent variable of the function may be another function of the C same group foregoing in the input data. For actual mapping of points C it is necessary to call the subroutine VAL2, which also returns the C first and second partial derivatives. Subroutine VAL1 may be called C several times. The groups in the class are indexed successively, C following the groups of the class defined during the previous C invocations. C C Input: C LUN... Logical unit number of the external input device C containing the input data. C ICLASS..Index of the class of the functions to be specified. The C classes are indexed by integers starting from 1. C NGROUP..Number of groups of functions to be specified during the C current invocation of VAL1. The groups of each class are C indexed by integers starting from 1. If some groups of C functions of the ICLASS-th class were specified in the C previous invocation of VAL1, the groups of functions now C read in are appended to them and are indexed following C them. C NFUNCT..Maximum number of functions to be specified for each C group. The actual number of specified functions may be C different for different groups. However, it must be less C than or equal to NFUNCT. C TFUNCT..Strings identifying the functions specified in the input C data. The function identified in the input data by string C TFUNCT(I) is associated with integer I. This integer C identifies what the function describes. C None of the input parameters are altered. C C No output. C C Common block: INCLUDE 'val.inc' C val.inc C All the storage locations of the common block are defined in this C subroutine. C C Subroutines and external functions required: * EXTERNAL CURVN1 EXTERNAL CURVB1 EXTERNAL VALB,SURFB1,VAL3B1,SORTV,READV C VALB... Block data subroutine of this file. C CURVN1 or CURVB1 (alternatives), SURFB1, VAL3B1, SNHCSH, VGEN, C TERMS, TRIDEC, TRISOL... Subroutine package 'FITPACK' C (file 'fit.for'). C SORTV, READV... This file. C C Date: 1999, September 15 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C LOGICAL WHAT INTEGER MCLASS,MGROUP,MFUNCT,LADR,MADR,MAXADR INTEGER KCLASS,KGROUP,KFUNCT,KADR,NVAR CHARACTER*3 TEXT REAL GROUP,SIGMA,POWERW INTEGER LX(3),LX1,LX2,LX3 EQUIVALENCE (LX(1),LX1),(LX(2),LX2),(LX(3),LX3) INTEGER NX(3),NX1,NX2,NX3 EQUIVALENCE (NX(1),NX1),(NX(2),NX2),(NX(3),NX3) INTEGER JADR(7),JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7 EQUIVALENCE (JADR(1),JADR1),(JADR(2),JADR2),(JADR(3),JADR3) EQUIVALENCE (JADR(4),JADR4),(JADR(5),JADR5),(JADR(6),JADR6) EQUIVALENCE (JADR(7),JADR7) INTEGER IGROUP,IFUNCT,IERR,I,J,L,N C C WHAT... Flag if the physical meaning of the functions is included C in the input data. C MCLASS,MGROUP,MFUNCT,LADR,MADR,MAXADR... Positions in the memory. C KCLASS,KGROUP,KFUNCT,KADR... Shifts in the memory. C NVAR... Number of the independent variables A1, A2, A3 of the C interpolated function W. C TEXT... String identifying the current group or the current C function. C GROUP...Index of the current group or power of the physical C quantity. C SIGMA...Tension factor. C POWERW..Given grid values (7) are raised to the (1/POWERW)-th C power immediately after reading and then interpolated. C LX=(LX1,LX2,LX3)... Addresses of auxiliary storage locations for C reordering the grid coordinates. C NX=(NX1,NX2,NX3)... Numbers of grid lines. C JADR=(JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7)... Addresses of C parameters describing the interpolated function (grid C coordinates, B-spline coefficients, B-spline basis C functions). C IGROUP,IFUNCT... Do loop variables. C IERR... Local variable to check the proper function of called C subroutines. C I,J,L,N... Local auxiliary variables. C C....................................................................... C C Flag if the physical meaning of the functions is included in the C input data: WHAT=.FALSE. DO 10 I=1,NFUNCT IF(TFUNCT(I).NE.' ') WHAT=.TRUE. 10 CONTINUE C C Positions in the memory: IF(ICLASS.GT.IPAR(0)) THEN C New class: MCLASS=IPAR(0) KCLASS=ICLASS-MCLASS ELSE C Old class: MCLASS=ICLASS KCLASS=0 END IF MGROUP=IPAR(MCLASS) KGROUP=KCLASS+NGROUP MFUNCT=IPAR(MGROUP) KFUNCT=KGROUP+NGROUP*NFUNCT MADR=IPAR(MFUNCT) KADR=NPAR-IPAR(IPAR(IPAR(IPAR(0)))) IF(KADR.LT.KFUNCT) GO TO 99 C Upper bound of the available memory MAXADR=MADR+KADR C C Movement in the memory: DO 11 I=IPAR(IPAR(IPAR(IPAR(0)))),MADR+1,-1 IPAR(I+KADR)=IPAR(I) 11 CONTINUE DO 12 I=MADR,IPAR(IPAR(IPAR(0)))+1,-1 IPAR(I+KFUNCT)=IPAR(I) 12 CONTINUE DO 13 I=IPAR(IPAR(IPAR(0))),MFUNCT+1,-1 IPAR(I+KFUNCT)=IPAR(I)+KADR 13 CONTINUE MADR=MADR+KFUNCT DO 14 I=MFUNCT,MGROUP+1,-1 IPAR(I+KGROUP)=IPAR(I)+KFUNCT 14 CONTINUE MFUNCT=MFUNCT+KGROUP DO 15 I=MGROUP,MCLASS+1,-1 IPAR(I+KCLASS)=IPAR(I)+KGROUP 15 CONTINUE MGROUP=MGROUP+KCLASS DO 16 I=0,MCLASS IPAR(I)=IPAR(I)+KCLASS 16 CONTINUE C New classes: DO 17 I=MCLASS+1,ICLASS IPAR(I)=IPAR(MCLASS) 17 CONTINUE C Number of previously stored groups of functions of the class IGROUP=IPAR(ICLASS)-IPAR(ICLASS-1) IPAR(ICLASS)=IPAR(ICLASS)+NGROUP C New groups: DO 18 I=MGROUP+1,MGROUP+NGROUP IPAR(I)=IPAR(I-1)+NFUNCT 18 CONTINUE C C Loop for groups of functions: GROUP=1. READ(LUN,*) TEXT,GROUP DO 90 IGROUP=IGROUP+1,IGROUP+NGROUP IF(TEXT.EQ.'END') THEN C 351 CALL ERROR('351 in VAL1: End of input functions encountered') C End of input functions encountered before all NGROUP C groups of functions are defined in the input data. ELSE IF(INT(GROUP+0.5).NE.IGROUP) THEN C 352 CALL ERROR('352 in VAL1: Improper index of the group') C Improper index of the group of input functions in the C input data. END IF C Loop for functions of the current group: DO 80 IFUNCT=1,NFUNCT C Physical meaning of the function: IF(WHAT) THEN GROUP=1. READ(LUN,*) TEXT,GROUP DO 21 I=1,NFUNCT IF(TFUNCT(I).EQ.TEXT) THEN MADR=MADR+2 IF(MADR.GT.MAXADR) GO TO 99 IPAR(MADR-1)=I RPAR(MADR)=GROUP GO TO 22 END IF 21 CONTINUE GO TO 81 22 CONTINUE ELSE MADR=MADR+2 IF(MADR.GT.MAXADR) GO TO 99 IPAR(MADR-1)=IFUNCT RPAR(MADR)=1. END IF C C Form of the function: LADR=MADR+1 MADR=MADR+4 IF(MADR.GT.MAXADR) GO TO 99 IPAR(LADR) =0 IPAR(LADR+1)=0 IPAR(LADR+2)=0 RPAR(MADR) =0. POWERW =1. READ(LUN,*) * IPAR(LADR),IPAR(LADR+1),IPAR(LADR+2),RPAR(MADR),POWERW SIGMA=RPAR(MADR) IF(POWERW.EQ.0.) THEN C 359 CALL ERROR('359 in VAL1: Zero POWERW in input data') C Zero POWERW in input data describing B-spline interpolation C of functions used to specify surfaces or material parameters C in the model or initial line or surface source. C See input data for one function, C line (2). END IF C Number of independent variables: NVAR=0 DO 23 I=LADR,LADR+2 IF(IPAR(I).GT.0) NVAR=NVAR+1 23 CONTINUE C C Numbers of grid coordinates: LADR=MADR+1 MADR=MADR+NVAR IF(MADR.GT.MAXADR) GO TO 99 IF(LADR.LE.MADR) THEN READ(LUN,*) (IPAR(I),I=LADR,MADR) END IF C C Reading grid coordinates: L=MAXADR+1 NVAR=0 DO 24 J=LADR,MADR N=IABS(IPAR(J)) IF(N.GT.0) THEN LADR=MADR+1 MADR=MADR+N IF(N.EQ.1) THEN IF(MADR.GE.L-1) GO TO 99 READ(LUN,*) RPAR(LADR) ELSE L=L-N IF(MADR+N.GE.L-1) GO TO 99 NVAR=NVAR+1 NX(NVAR)=N LX(NVAR)=L JADR(NVAR)=LADR READ(LUN,*) (RPAR(I),I=MADR+1,MADR+N) CALL SORTV(N,RPAR(MADR+1),RPAR(LADR),IPAR(L)) END IF END IF 24 CONTINUE DO 25 I=NVAR+1,3 NX(I)=1 LX(I)=L-1 IPAR(L-1)=1 25 CONTINUE C C Reading grid values: JADR4=MADR+1 MADR=MADR+NX1*NX2*NX3 IF(MADR.GE.L) GO TO 99 CALL READV(LUN,NX1,NX2,NX3,IPAR(LX1),IPAR(LX2),IPAR(LX3), * RPAR(JADR4),POWERW) C C Computing B-spline under tension expansion coefficients: IF(NVAR.LE.0) THEN C No independent variable: CONTINUE ELSE C Size of the temporary storage location N=3*MAX0(NX1,NX2,NX3) JADR5=MADR+1 MADR=MADR+5*NX1 IF(MADR+N.GT.MAXADR) GO TO 99 C IERR enables to check for the proper function of subroutines C called IERR=1 IF(NVAR.EQ.1) THEN C One independent variable: C Two alternatives: Hermite or B-spline representations C may be used for the 1-D interpolation. Just one of the C following two statements must be supplied by '*' in the C first column: C First statement - Hermite representation: * CALL CURVN1(NX1,RPAR(JADR1),RPAR(JADR4), * * RPAR(JADR5),RPAR(MADR+1),SIGMA,IERR) C Second statement - B-spline representation: CALL CURVB1(NX1,RPAR(JADR1),RPAR(JADR4),RPAR(JADR4), * RPAR(JADR5),RPAR(MADR+1),SIGMA,IERR) C Do not forget to supply '*' into the first column of the C corresponding statement in subroutine VAL2. ELSE JADR6=MADR+1 MADR=MADR+5*NX2 IF(MADR+N.GT.MAXADR) GO TO 99 IF(NVAR.EQ.2) THEN C Two independent variables: CALL SURFB1(NX1,NX2,RPAR(JADR1),RPAR(JADR2), * RPAR(JADR4),NX1,RPAR(JADR4), * RPAR(JADR5),RPAR(JADR6),RPAR(MADR+1),SIGMA,IERR) ELSE C Three independent variables: JADR7=MADR+1 MADR=MADR+5*NX3 IF(MADR+N.GT.MAXADR) GO TO 99 CALL VAL3B1(NX1,NX2,NX3,RPAR(JADR1),RPAR(JADR2), * RPAR(JADR3),RPAR(JADR4),NX1,NX2,RPAR(JADR4), * RPAR(JADR5),RPAR(JADR6),RPAR(JADR7), * RPAR(MADR+1),SIGMA,IERR) END IF END IF IF(IERR.NE.0) THEN C 353 CALL ERROR('353 in VAL1: Strange error') C This error in the input functions should not appear. C Contact the authors. END IF END IF C Coefficients are evaluated C MFUNCT=MFUNCT+1 IPAR(MFUNCT)=MADR 80 CONTINUE GROUP=1. READ(LUN,*) TEXT,GROUP C End of loop for functions C C The remaining functions of the current group are not defined by C the input data: 81 CONTINUE DO 82 I=IFUNCT,NFUNCT MFUNCT=MFUNCT+1 IPAR(MFUNCT)=MADR 82 CONTINUE 90 CONTINUE C End of loop for groups of functions C IF(TEXT.NE.'END') THEN C 354 CALL ERROR('354 in VAL1: Input functions not properly ended') C Read in input data describing functions are not properly ended. END IF C C Movement in the memory: KADR=MAXADR-MADR DO 91 I=MAXADR+1,NPAR IPAR(I-KADR)=IPAR(I) 91 CONTINUE DO 92 I=MFUNCT+1,IPAR(IPAR(IPAR(0))) IPAR(I)=IPAR(I)-KADR 92 CONTINUE RETURN C 99 CONTINUE C 355 CALL ERROR('355 in VAL1: Insufficient memory in /VALC/') C Insufficient memory for the input data in common block /VALC/. C The dimension NPAR of array IPAR (or RPAR) must be enlarged. C See the block data subroutine VALB. END C C----------------------------------------------------------------------- C C C SUBROUTINE SORTV(NX,X1,X2,IX) INTEGER NX,IX(NX) REAL X1(NX),X2(NX) C C This subroutine is an auxiliary routine to VAL1. It reorders the C input grid coordinates to be ascending. C C Auxiliary storage locations INTEGER I,J C DO 3 J=1,NX IX(J)=1 DO 1 I=1,J-1 IF(X1(J).EQ.X1(I)) GO TO 9 IF(X1(J).GT.X1(I)) IX(J)=IX(J)+1 1 CONTINUE DO 2 I=J+1,NX IF(X1(J).EQ.X1(I)) GO TO 9 IF(X1(J).GT.X1(I)) IX(J)=IX(J)+1 2 CONTINUE 3 CONTINUE DO 4 J=1,NX X2(IX(J))=X1(J) 4 CONTINUE RETURN C 9 CONTINUE C 356 CALL ERROR('356 in SORTV in VAL1: Identical grid coordinates') C Two identical grid coordinates encountered in the input data. END C C----------------------------------------------------------------------- C C C SUBROUTINE READV(LUN,NX1,NX2,NX3,IX1,IX2,IX3,VAL,POWERW) INTEGER LUN,NX1,NX2,NX3,IX1(NX1),IX2(NX2),IX3(NX3) REAL VAL(NX1,NX2,NX3),POWERW C C This subroutine is an auxiliary routine to VAL1. It reads from the C input data the values given at grid points. C C Auxiliary storage locations INTEGER I1,I2,I3 REAL AUX1 C READ(LUN,*) (((VAL(IX1(I1),IX2(I2),IX3(I3)),I1=1,NX1), * I2=1,NX2),I3=1,NX3) IF(POWERW.NE.1.) THEN AUX1=1./POWERW DO 3 I3=1,NX3 DO 2 I2=1,NX2 DO 1 I1=1,NX1 VAL(IX1(I1),IX2(I2),IX3(I3))= * VAL(IX1(I1),IX2(I2),IX3(I3))**AUX1 1 CONTINUE 2 CONTINUE 3 CONTINUE END IF C RETURN END C C======================================================================= C C C SUBROUTINE VAL2(ICLASS,IGROUP,NFUNCT,COOR,F,POWER) INTEGER ICLASS,IGROUP,NFUNCT REAL COOR(3),F(10,NFUNCT),POWER(NFUNCT) C C This subroutine evaluates the function value, the three first partial C derivatives and the six second partial derivatives of a given function C at a given point. C C Input: C ICLASS..Index of the class of the required functions. The classes C are indexed by integers starting from 1. C IGROUP..Index of the group of the required functions. The groups C of each class are indexed by integers starting from 1. C NFUNCT..Number of the required functions. All functions belonging C to the IGROUP-th group of the ICLASS-th class and defined C by the input data must be required. The functions defined C by the input data (see subroutine VAL1) are one-to-one C corresponding to the integers which identify what the C function describes. The position of each evaluated C function in the output array F (see below) is determined C by this integer. That is why NFUNCT must be greater than C or equal to the greatest of these integers. The required C functions not defined by the input data are defined on the C output of this subroutine and are zero. C COOR... Array containing coordinates X1, X2, X3 of the given point C None of the input parameters are altered. C C Output: C F... Array containing, in each its column, function value, the C first and second partial derivatives of the corresponding C evaluated function in the order F, F1, F2, F3, F11, F12, C F22, F13, F23, F33. C POWER...The specified function is equal to the POWER-th power of C the corresponding physical quantity. The zero value of C the POWER indicates that the corresponding function is not C defined by the input data. C C Common block: INCLUDE 'val.inc' C val.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: * EXTERNAL CURV2D EXTERNAL CURVBD EXTERNAL SURFBD,VAL3BD C CURV2D or CURVBD (alternatives), SURFBD, VAL3BD, SNHCSH, DSPLNZ, C INTRVL... Subroutine package 'FITPACK' (file 'fit.for'). C C Date: 1995, March 28 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C The evaluated function has the form of C F(X1,X2,X3) = W(A1,A2,A3) - B1 - B2 - B3 , C C Its first derivatives are C dF dW dAk dB1 dB2 dB3 C --- = --- * --- - --- - --- - --- , C dXi dAk dXi dXi dXi dXi C C Its second derivatives are C d2 F d W d2 Ak d2 W dAk dAj d2 B1 d2 B3 C ------- = ---*------- + -------*---*--- - ------- - ... - -------. C dXi dXm dAk dXi dXm dAk dAj dXi dXm dXi dXm dXi dXm C C....................................................................... C INTEGER JGROUP,LFUNCT,MFUNCT,JFUNCT,LADR,MADR,IADR,IVAL INTEGER NVAR,IVAR(3),JVAR,KVAR INTEGER NX(3),NX1,NX2,NX3 EQUIVALENCE (NX(1),NX1),(NX(2),NX2),(NX(3),NX3) REAL XX(3),XX1,XX2,XX3 EQUIVALENCE (XX(1),XX1),(XX(2),XX2),(XX(3),XX3) CV3 REAL R1,R2,R3 INTEGER JADR(7),JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7 EQUIVALENCE (JADR(1),JADR1),(JADR(2),JADR2),(JADR(3),JADR3) EQUIVALENCE (JADR(4),JADR4),(JADR(5),JADR5),(JADR(6),JADR6) EQUIVALENCE (JADR(7),JADR7) REAL SIGMA,W(10),AUX1,AUX2 INTEGER I,J,K,M,N,ISYM(3,3) DATA ISYM/5,6,8,6,7,9,8,9,10/ C C JGROUP..Address of the IGROUP-th group of the ICLASS-th class. C LFUNCT,MFUNCT,JFUNCT... Addresses of the first, last and arbitrary C functions of the group. C LADR,MADR,IADR... Addresses of the first, last and arbitrary C parameters of the current function. C IVAL... Index of the function F being currently evaluated. C NVAR,IVAR(3),JVAR,KVAR... Number and types of the independent C variables A1, A2, A3 of the interpolated function W. C NX=(NX1,NX2,NX3)... Numbers of grid lines. C XX=(XX1,XX2,XX3),R1,R2,R3... Values of independent variables A1, C A2, A3 of function W. C JADR=(JADR1,JADR2,JADR3,JADR4,JADR5,JADR6,JADR7)... Addresses of C parameters describing the interpolated function (grid C coordinates, B-spline coefficients, B-spline basis C functions). C SIGMA...Tension factor. C W... Array for the value, the first and second partial C derivatives of function W. C AUX1,AUX2,I,J,K,M,N... Local auxiliary variables. C ISYM... Storage of the symmetric 3*3 matrix. C C....................................................................... C C The default value of the function is the zero function. C Loop for the functions to be evaluated: DO 12 J=1,NFUNCT DO 11 I=1,10 F(I,J)=0. 11 CONTINUE POWER(J)=0. 12 CONTINUE *V CALL VAR1() C IF(ICLASS.LT.1.OR.IPAR(0).LT.ICLASS) THEN C 357 WRITE(*,'(2(A,I10))') ' CLASS=',ICLASS,', GROUP =',IGROUP CALL ERROR('357 in VAL2: Incorrect index of the class') C The index of the class of the functions to be evaluated is zero, C negative or greater than the number of classes defined. END IF JGROUP=IPAR(ICLASS-1)+IGROUP IF(IGROUP.LT.1.OR.IPAR(ICLASS).LT.JGROUP) THEN C 358 WRITE(*,'(2(A,I10))') ' CLASS=',ICLASS,', GROUP =',IGROUP CALL ERROR('358 in VAL2: Incorrect index of the group') C The index of the group of the functions to be evaluated is zero, C negative or greater than the number of groups defined within the C given class. END IF LFUNCT=IPAR(JGROUP-1)+1 MFUNCT=IPAR(JGROUP) MADR =IPAR(LFUNCT-1) C C Loop for functions F being evaluated: DO 90 JFUNCT=LFUNCT,MFUNCT C Starting and end addresses of the parameters describing the C function LADR=MADR+1 MADR=IPAR(JFUNCT) IF(LADR.LE.MADR) THEN C Index of function F being currently evaluated IVAL=IPAR(LADR) C Power of the corresponding physical quantity POWER(IVAL)=RPAR(LADR+1) C Tension factor SIGMA=RPAR(LADR+5) C C The number, types and values of the independent variables Ai C of function W being interpolated, and the functions Bi being C subtracted from the evaluated function: C Initial address IADR=LADR+6 C Initial number of the independent variables NVAR=0 JADR1=0 JADR2=0 JADR3=0 JADR4=0 C Loop for the possible independent variables: DO 20 M=LADR+2,LADR+4 C Type of the possible independent variable: J=IPAR(M) IF(J.NE.0) THEN IF(J.GT.0) THEN N=IABS(IPAR(IADR)) IF(N.GE.2) THEN NVAR=NVAR+1 NX(NVAR)=N IF(J.LE.3) THEN IVAR(NVAR)=J XX(NVAR)=COOR(J) ELSE K=IPAR(IPAR(LFUNCT+J-5)+1) IVAR(NVAR)=K+3 XX(NVAR)=F(1,K) END IF ELSE IF(N.EQ.1) THEN JADR(NVAR+1)=JADR(NVAR+1)+1 END IF IADR=IADR+1 ELSE C Subtracting certain functions from function F being C evaluated: IF(J.GE.-3) THEN C Subtracting a coordinate: F(1,IVAL)=F(1,IVAL)-COOR(-J) F(1-J,IVAL)=F(1-J,IVAL)-1. ELSE C Subtracting another function F: K=IPAR(IPAR(LFUNCT-J-5)) DO 19 I=1,10 F(I,IVAL)=F(I,IVAL)-F(I,K) 19 CONTINUE *V CALL VAR4(0,-1.) *V CALL VAR5(IVAL,K) END IF END IF END IF 20 CONTINUE C CV3 Lines denoted by 'CV3' in the first 3 columns are related to CV3 the kind of extrapolation outside the grid used in version 3 CV3 (January 1991) and older. In those versions, the first and CV3 second derivatives were incorrect outside the grid. CV3 If removing 'CV3' and 'CV3-V' from the executable statements, CV3 the kind of extrapolation from ver.3 is restored. Then, the CV3 first derivatives are correctly evaluated (unlike in ver.3), CV3 but the second derivatives are incorrect (as in ver.3). CV3 Similarly, variations of functional values are correct, and CV3 first variations of first derivatives are incorrect. CV3 C Interpolation of function W: JADR1=IADR+JADR1 CV3-V CALL VAR4(0,1.) IF(NVAR.LE.0) THEN C No independent variable: W(1)=RPAR(JADR1) *V CALL VAR2(1,1.,0.,0.,0.) *V CALL VAR3(JADR1-1) ELSE JADR2=JADR1+NX1+JADR2 CV3 R1=XX1 CV3 IF(XX1.LT.RPAR(JADR1)) THEN CV3 XX1=RPAR(JADR1) CV3 ELSE IF(XX1.GT.RPAR(JADR2-1)) THEN CV3 XX1=RPAR(JADR2-1) CV3 END IF CV3 R1=R1-XX1 IF(NVAR.EQ.1) THEN C One independent variable: JADR3=JADR2+NX1 C Two alternatives: Hermite or B-spline representations C may be used for the 1-D interpolation. Just one of the C following two statements must be supplied by '*' in the C first column: C First statement - Hermite representation: * CALL CURV2D(XX1,W(1),W(2),W(5),NX1, * * RPAR(JADR1),RPAR(JADR2),RPAR(JADR3),SIGMA) C Second statement - B-spline representation: CALL CURVBD(XX1,W(1),W(2),W(5),NX1, * RPAR(JADR1),RPAR(JADR2),RPAR(JADR3),SIGMA) C Do not forget to supply '*' into the first column of the C corresponding statement in subroutine VAL1. *V CALL VAR3(JADR2-1) ELSE JADR3=JADR2+NX2+JADR3 CV3 R2=XX2 CV3 IF(XX2.LT.RPAR(JADR2)) THEN CV3 XX2=RPAR(JADR2) CV3 ELSE IF(XX2.GT.RPAR(JADR3-1)) THEN CV3 XX2=RPAR(JADR3-1) CV3 END IF CV3 R2=R2-XX2 IF(NVAR.EQ.2) THEN C Two independent variables: JADR4=JADR3+NX1*NX2 JADR5=JADR4+5*NX1 CALL SURFBD(XX1,XX2,W(1),W(2),W(3),W(5),W(6),W(7), * NX1,NX2,RPAR(JADR1),RPAR(JADR2),RPAR(JADR3), * RPAR(JADR4),RPAR(JADR5),SIGMA) *V CALL VAR3(JADR3-1) ELSE C Three independent variables: JADR4=JADR3+NX3+JADR4 JADR5=JADR4+NX1*NX2*NX3 JADR6=JADR5+5*NX1 JADR7=JADR6+5*NX2 CV3 R3=XX3 CV3 IF(XX3.LT.RPAR(JADR3)) THEN CV3 XX3=RPAR(JADR3) CV3 ELSE IF(XX3.GT.RPAR(JADR4-1)) THEN CV3 XX3=RPAR(JADR4-1) CV3 END IF CV3 R3=R3-XX3 CALL VAL3BD(XX1,XX2,XX3,W(1),W(2),W(3),W(4),W(5),W(6), * W(7),W(9),W(10),W(8),NX1,NX2,NX3, * RPAR(JADR1),RPAR(JADR2),RPAR(JADR3),RPAR(JADR4), * RPAR(JADR5),RPAR(JADR6),RPAR(JADR7),SIGMA) *V CALL VAR3(JADR4-1) CV3 W(1)=W(1)+W(4)*R3 CV3 IF(R1.EQ.0.) W(2)=W(2)+W(8)*R3 CV3 IF(R2.EQ.0.) W(3)=W(3)+W(9)*R3 CV3 IF(R3.EQ.0.) W(4)=W(4)+W(8)*R1+W(9)*R2 CV3-V CALL VAR4(13,R3) END IF CV3 W(1)=W(1)+W(3)*R2 CV3 IF(R1.EQ.0.) W(2)=W(2)+W(6)*R2 CV3 IF(R2.EQ.0.) W(3)=W(3)+W(6)*R1 CV3-V CALL VAR4(9,R2) END IF CV3 W(1)=W(1)+W(2)*R1 CV3-V CALL VAR4(5,R1) END IF CV3-V CALL VAR5(0,0) C Function W is evaluated C C Evaluation of function f: C Functional value (zero derivative) F(1,IVAL)=F(1,IVAL)+W(1) *V CALL VAR4(0,0.) *V CALL VAR4(1,1.) C Loop for the summation index K: DO 39 K=1,NVAR KVAR=IVAR(K) IF(KVAR.LE.3) THEN C First derivatives - first term on R.H.S. F(1+KVAR,IVAL)=F(1+KVAR,IVAL)+W(1+K) C Second derivatives - second term on R.H.S. (the first term C vanishes in this case) - loop for the summation index J: DO 32 J=1,NVAR JVAR=IVAR(J) IF(JVAR.LE.3) THEN IF(JVAR.LE.KVAR) THEN N=ISYM(JVAR,KVAR) F(N,IVAL)=F(N,IVAL)+W(ISYM(J,K)) END IF ELSE JVAR=JVAR-3 AUX1=W(ISYM(J,K)) DO 31 I=1,JVAR N=ISYM(I,JVAR) F(N,IVAL)=F(N,IVAL)+AUX1*F(1+I,JVAR) 31 CONTINUE END IF 32 CONTINUE *V CALL VAR4(4*K+1+KVAR,1.) ELSE KVAR=KVAR-3 DO 33 I=2,4 *V CALL VAR4(4*K+I,F(I,KVAR)) 33 CONTINUE END IF 39 CONTINUE *V CALL VAR5(IVAL,0) C Loop for the summation index K: DO 49 K=1,NVAR KVAR=IVAR(K) IF(KVAR.GT.3) THEN KVAR=KVAR-3 *V CALL VAR4(0,W(1+K)) C First and second derivatives - first terms on R.H.S. DO 44 I=2,10 F(I,IVAL)=F(I,IVAL)+W(1+K)*F(I,KVAR) 44 CONTINUE C Second derivatives - second term on R.H.S. - C loop for the summation index J: DO 48 J=1,NVAR JVAR=IVAR(J) IF(JVAR.LE.3) THEN AUX1=W(ISYM(J,K)) DO 45 I=1,KVAR N=ISYM(I,KVAR) F(N,IVAL)=F(N,IVAL)+AUX1*F(1+I,KVAR) 45 CONTINUE ELSE JVAR=JVAR-3 AUX1=W(ISYM(J,K)) DO 47 M=1,3 AUX2=AUX1*F(1+M,JVAR) DO 46 I=1,M N=ISYM(I,M) F(N,IVAL)=F(N,IVAL)+AUX2*F(1+I,KVAR) 46 CONTINUE *V CALL VAR4(1+M,AUX2) 47 CONTINUE END IF 48 CONTINUE *V CALL VAR5(IVAL,KVAR) END IF 49 CONTINUE C END IF 90 CONTINUE C End of loop for evaluated functions F C RETURN END C C======================================================================= C