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 By Vlastislav Cerveny, Ludek Klimes, Ivan Psencik
C
C This file consists of the following subroutines:
C VALB... Block data subroutine defining common block /VALC/ to
C store the data describing the interpolated functions.
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 SORTV...Auxiliary subroutine to subroutine VAL1.
C READV...Auxiliary subroutine to subroutine VAL1.
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 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 Attention:
C The lines denoted by '*V' in the first 2 columns should be removed
C from the subroutine VAL2 before compilation.
C If the model variations with respect to the model parameters were
C to be calculated, these lines would remain while replacing each
C '*V' in the first columns by spaces. In such a case, subroutines
C VAR1, VAR2, VAR3, VAR4 and VAR5 of the file 'var.for' would be
C called to handle the variations.
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 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', 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 apriori 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 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 (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 Storage in the memory:
C The parameters describing the interpolated functions are stored
C in common block /VALC/ defined in the following subroutine:
C ------------------------------------------------------------------
BLOCK DATA VALB
INCLUDE 'val.inc'
DATA IPAR(0)/0/
END
C ------------------------------------------------------------------
C
C Date: 1996, September 30
C Coded by Ludek Klimes
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 succesively,
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 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 Error messages:
C 351... End of input functions encountered:
C End of input functions encountered before all ngroup
C groups of functions are defined in the input data.
C 352... Improper index of the group:
C Improper index of the group of input functions in the
C input data.
C 353... Strange error:
C This error in the input functions should not appear.
C Contact the authors.
C 354... Input functions not properly ended:
C Read in input data describing functions are not properly
C ended.
C 355... Insufficient memory in /VALC/:
C Insufficient memory for the input data in common block
C /VALC/. The dimension npar of array IPAR (or RPAR) must
C be enlarged. See the block data subroutine VALB.
C 356... Identical grid coordinates:
C Two identical grid coordinates encountered in the input
C data.
C
C Date: 1996, July 22
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:
READ(LUN,*) TEXT,GROUP
DO 90 IGROUP=IGROUP+1,IGROUP+NGROUP
IF(TEXT.EQ.'END') THEN
PAUSE 'Error 351 in VAL1: End of input functions encountered'
STOP
ELSE IF(INT(GROUP+0.5).NE.IGROUP) THEN
PAUSE 'Error 352 in VAL1: Improper index of the group'
STOP
END IF
C Loop for functions of the current group:
DO 80 IFUNCT=1,NFUNCT
C Physical meaning of the function:
IF(WHAT) THEN
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)
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
PAUSE 'Error 353 in VAL1: Strange error'
STOP
END IF
END IF
C Coefficients are evaluated
C
MFUNCT=MFUNCT+1
IPAR(MFUNCT)=MADR
80 CONTINUE
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
PAUSE 'Error 354 in VAL1: Input functions not properly ended'
STOP
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
PAUSE 'Error 355 in VAL1: Insufficient memory in /VALC/'
STOP
END
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
PAUSE 'Error 356 in SORTV in VAL1: Identical grid coordinates'
STOP
END
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
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 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 Error messages:
C 357... Incorrect index of the class:
C The index of the class of the functions to be evaluated is
C zero, negative or greater than the number of classes
C defined.
C 358... Incorrect index of the group:
C The index of the group of the functions to be evaluated is
C zero, negative or greater than the number of groups
C defined within the given class.
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
CALL VAR1()
C
IF(ICLASS.LT.1.OR.IPAR(0).LT.ICLASS) THEN
WRITE(*,'(2(A,I10))') ' CLASS=',ICLASS,', GROUP =',IGROUP
PAUSE 'Error 357 in VAL2: Incorrect index of the class'
STOP
END IF
JGROUP=IPAR(ICLASS-1)+IGROUP
IF(IGROUP.LT.1.OR.IPAR(ICLASS).LT.JGROUP) THEN
WRITE(*,'(2(A,I10))') ' CLASS=',ICLASS,', GROUP =',IGROUP
PAUSE 'Error 358 in VAL2: Incorrect index of the group'
STOP
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
CALL VAR4(0,-1.)
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)
CALL VAR2(1,1.,0.,0.,0.)
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.
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)
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)
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)
CALL VAR4(0,0.)
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
CALL VAR4(4*K+1+KVAR,1.)
ELSE
KVAR=KVAR-3
DO 33 I=2,4
CALL VAR4(4*K+I,F(I,KVAR))
33 CONTINUE
END IF
39 CONTINUE
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
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
CALL VAR4(1+M,AUX2)
47 CONTINUE
END IF
48 CONTINUE
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