C 
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