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 COLUMN SHOULD BE REMOVED
C     FROM THE SUBROUTINE VAL 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 COLUMN BY A SPACE.  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
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 (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 FUNCTION W AT GRID POINTS. FUNCTION VALUE W(I,J,K)
C     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
      INTEGER NPAR
      PARAMETER (NPAR=8191)
      INTEGER IPAR(0:NPAR)
      REAL RPAR(0:NPAR)
      EQUIVALENCE (IPAR,RPAR)
      COMMON/VALC/IPAR
      SAVE /VALC/
      DATA IPAR(0)/0/
      END
C     ------------------------------------------------------------------
C     COMMON BLOCK /VALC/ IS INCLUDED IN SUBROUTINES VAL1 AND VAL2.
C     THE PARAMETERS ARE STORED SEQUENTIALLY IN ONE ARRAY REGARDLESS OF
C     THE FACT WHETHER THEY ARE OF TYPE INTEGER OR REAL.  THE INDIVIDUAL
C     NUMERIC STORAGE UNITS OF THE ARRAY ARE INDEXED STARTING FROM 0 AND
C     ARE NAMED IPAR(I) OR RPAR(I) DEPENDING ON THE TYPE OF A PARAMETER.
C     THE INDEX OF THE LAST ALLOCATED NUMERIC STORAGE UNIT IS NPAR (SEE
C     (F) BELOW).  IF NPAR IS CHANGED (SEE THE THIRD STATEMENT OF THE
C     ABOVE BLOCK DATA SUBROUTINE VALB), IT MUST BE ADJUSTED IN
C     SUBROUTINES VAL1 AND VAL2, TOO.
C
C     COMMON BLOCK /VALC/ CAN BE DIVIDED INTO TWO PARTS.  THE FIRST PART
C     (SEE (A)-(D) BELOW) OF THE COMMON BLOCK /VALC/ CONTAINS INTEGERS.
C     THE NUMBER OF THESE INTEGERS IN THE FIRST PART EQUALS TO
C     1+MCLASS+NG+NF, WHERE MCLASS IS THE TOTAL NUMBER OF CLASSES, NG IS
C     THE TOTAL NUMBER OF ALL GROUPS OF ALL CLASSES, NF IS THE TOTAL
C     NUMBER OF ALL FUNCTIONS OF ALL GROUPS OF ALL CLASSES. THE FIRST
C     PART OF THE COMMON BLOCK SPECIFIES THE DIVISION OF THE COMMON
C     BLOCK /VALC/ INTO PARAMETERS DESCRIBING INDIVIDUAL FUNCTIONS.
C     THE SECOND PART (SEE (E) BELOW) OF THE COMMON BLOCK /VALC/
C     CONTAINS THE PARAMETERS DESCRIBING INDIVIDUAL FUNCTIONS, STORED
C     SUCCESSSIVELY FOR THE FIRST, SECOND, ..., LAST FUNCTION OF THE
C     FIRST GROUP OF THE FIRST CLASS, FOR THE FIRST, SECOND, ..., LAST
C     FUNCTION OF THE SECOND GROUP OF THE FIRST CLASS, ..., FOR THE
C     FIRST, SECOND, ..., LAST FUNCTION OF THE FIRST GROUP OF THE SECOND
C     CLASS, ..., FOR THE FIRST, SECOND, ..., LAST FUNCTION OF THE LAST
C     GROUP OF THE LAST CLASS.  THE SECOND PART OF THE COMMON BLOCK
C     /VALC/ CONTAINS BOTH INTEGER AND REAL PARAMETERS.
C     FOR AN EXAMPLE REFER TO THE SAMPLE INPUT DATA FOR THE MODEL.
C
C     (A) NUMERIC STORAGE UNIT IPAR(0) CONTAINS THE NUMBER MCLASS OF
C     CLASSES.
C
C     (B) CLASSES:
C     SUBSEQUENT NUMERIC STORAGE UNITS (IPAR(I),I=1,MCLASS) CORRESPOND
C     TO THE INDIVIDUAL CLASSES.  ONE NUMERIC STORAGE UNIT CORRESPONDS
C     TO ONE CLASS.  THE NUMERIC STORAGE UNIT CORRESPONDING TO A CLASS
C     CONTAINS THE INDEX OF THE NUMERIC STORAGE UNIT CORRESPONDING TO
C     THE LAST GROUP OF THE CLASS.
C
C     (C) GROUPS:
C     SUBSEQUENT NUMERIC STORAGE UNITS (IPAR(I),I=MCLASS+1,MCLASS+NG)
C     CORRESPOND TO THE INDIVIDUAL GROUPS.  THE TOTAL NUMBER OF THESE
C     STORAGE UNITS IS NG.  ONE NUMERIC STORAGE UNIT CORRESPONDS TO ONE
C     GROUP.  THE NUMERIC STORAGE UNIT CORRESPONDING TO A GROUP CONTAINS
C     THE INDEX OF THE NUMERIC STORAGE UNIT CORRESPONDING TO THE LAST
C     FUNCTION OF THE GROUP.
C
C     (D) FUNCTIONS:
C     SUBSEQUENT NUMERIC STORAGE UNITS (IPAR(I),I=MCLASS+NG+1,
C     MCLASS+NG+NF) CORRESPOND TO THE INDIVIDUAL FUNCTIONS.  THE TOTAL
C     NUMBER OF THESE STORAGE UNITS IS NF.  ONE NUMERIC STORAGE UNIT
C     CORRESPONDS TO ONE FUNCTION.  THE NUMERIC STORAGE UNIT
C     CORRESPONDING TO A FUNCTION CONTAINS THE INDEX OF THE NUMERIC
C     STORAGE UNIT CORRESPONDING TO THE LAST PARAMETER THAT DESCRIBES
C     THE FUNCTION.  THE FUNCTIONS ARE STORED IN THE ORDER IN WHICH THEY
C     HAVE BEEN SPECIFIED BY THE INPUT DATA.  SINCE THE NUMBER OF
C     FUNCTIONS IN THE INPUT DATA CORRESPODING TO ONE GROUP MAY BE LESS
C     THAN THEIR MAXIMUM NUMBER NFUNCT (INPUT PARAMETER OF SUBROUTINE
C     VAL1), SOME NUMERIC STORAGE UNITS RELEVANT TO FUNCTIONS MAY NOT
C     CORRESPOND TO THE FUNCTIONS SPECIFIED BY THE INPUT DATA.  EACH
C     SUCH NUMERIC STORAGE UNIT CONTAINS THE SAME ADDRESS AS THE
C     PREVIOUS NUMERIC STORAGE UNIT (I.E.  THE CORRESPONDING FUNCTION IS
C     SPECIFIED BY NO PARAMETER) AND HAS NO INFLUENCE ON THE FUNCTION
C     EVALUATION.
C
C     (E) PARAMETERS OF FUNCTIONS:
C     SUBSEQUENT NUMERIC STORAGE UNITS (IPAR(I),I=MCLASS+NG+NF+1,
C     MCLASS+NG+NF+NP) CONTAIN THE PARAMETERS DESCRIBING INDIVIDUAL
C     FUNCTIONS.  HERE WE HAVE DENOTED BY NP THE TOTAL NUMBER OF THESE
C     STORAGE UNITS.  ANY NUMBER OF NUMERIC STORAGE UNITS MAY CORRESPOND
C     TO THE PARAMETERS OF ONE FUNCTION.  THE FIRST NUMERIC STORAGE UNIT
C     CORRESPONDING TO A FUNCTION CONTAINS THE INTEGER THAT IDENTIFIES
C     THE PHYSICAL MEANING OF THE FUNCTION, E.G.  ITS VALUES MAY
C     IDENTIFY P-WAVE VELOCITY, S-WAVE VELOCITY, DENSITY, ETC.  THE
C     SECOND NUMERIC STORAGE UNIT CONTAINS THE POWER OF THE PHYSICAL
C     QUANTITY, SEE 'INPUT DATA FOR ONE FUNCTION' (1).  IF ITEM (1) OF
C     'INPUT DATA FOR ONE FUNCTION' IS OMITTED, THE TWO FIRST NUMERIC
C     STORAGE UNITS CORRESPONDING TO THE PARAMETERS OF THE FUNCTION
C     CONTAIN 1 (THE FIRST INTEGER, THE SECOND REAL).  THE SUBSEQUENT
C     NUMERIC STORAGE UNITS CONTAIN 'INPUT DATA FOR ONE FUNCTION' (2),
C     (3), (4), (5) AND (6), STORED IN THE SAME AMOUNT AS READ FROM THE
C     INPUT DATA.  DATA (4), (5), (6) AND (7) ARE REORDERED TO RENDER
C     THE GRID COORDINATES IN ASCENDING ORDER.  INSTEAD OF THE GRID
C     VALUES (7), THE COEFFICIENTS DESCRIBING THE FUNCTION IN TERMS OF A
C     B-SPLINE UNDER TENSION BASIS ARE STORED (ONLY IF SUBROUTINE CURVN1
C     (HERMITE REPRESENTATION) IS USED, THE GRID VALUES (7) ARE STORED).
C     THESE PARAMETERS ARE FOLLOWED BY ARRAYS OF LENGTH 5*NX(1), ...,
C     5*NX(NVAR) RESPECTIVELY, CONTAINING THE B-SPLINE UNDER TENSION
C     BASIS DATA COMPUTED FOR THE PROJECTIONS OF THE GRID ONTO THE AXES
C     OF INDIVIDUAL INDEPENDENT VARIABLES (ONLY IF SUBROUTINE CURVN1
C     (HERMITE REPRESENTATION) IS USED, THE FIRST NX(1) NUMERIC STORAGE
C     UNITS CONTAIN THE SECOND DERIVATIVES AND THE FOLLOWING 4*NX(1)
C     NUMERIC STORAGE UNITS ARE UNDEFINED.
C
C     (F) UNDEFINED PART OF THE COMMON BLOCK:
C     SUBSEQUENT NUMERIC STORAGE UNITS (IPAR(I),I=MCLASS+NG+NF+NP+1,
C     NPAR) ARE UNDEFINED.
C
C     MEMORY MODEL OF THE COMMON BLOCK /VALC/:
C     ------------------------------------------------------------------
C     ADDRESS:     0    JCLASS-1 JCLASS JGROUP-1 JGROUP JFUNCT-1 JFUNCT
C     VALUE:     MCLASS LGROUP-1 MGROUP LFUNCT-1 MFUNCT  LADR-1   MADR
C     ..................................................................
C     FOR THE MEANING OF INDIVIDUAL ITEMS IN THIS TABLE SEE BELOW.
C     ------------------------------------------------------------------
C     THE WAY OF ACCESS TO IFUNCT-TH FUNCTION OF IGROUP-TH GROUP OF THE
C     ICLASS-TH CLASS:
C     ADDRESS OF LAST CLASS..................... MCLASS=IPAR(0)
C     ADDRESS OF ICLASS-TH CLASS................ JCLASS=ICLASS
C             (IT MUST BE: 1.LE.JCLASS.LE.MCLASS)
C     ADDRESS OF FIRST GROUP OF THE CLASS....... LGROUP=IPAR(JCLASS-1)+1
C     ADDRESS OF LAST GROUP OF THE CLASS........ MGROUP=IPAR(JCLASS)
C     ADDRESS OF IGROUP-TH GROUP OF THE CLASS... JGROUP=LGROUP-1+IGROUP
C             (IT MUST BE: LGROUP.LE.JGROUP.LE.MGROUP)
C     ADDRESS OF FIRST FUNCTION OF THE GROUP.... LFUNCT=IPAR(JGROUP-1)+1
C     ADDRESS OF LAST FUNCTION OF THE GROUP..... MFUNCT=IPAR(JGROUP)
C     ADDRESS OF IFUNCT-TH FUNCTION OF THE GROUP JFUNCT=LFUNCT-1+IFUNCT
C             (IT MUST BE: LFUNCT.LE.JFUNCT.LE.MFUNCT)
C     ADDRESS OF FIRST FUNCTION PARAMETER....... LADR=IPAR(JFUNCT-1)+1
C     ADDRESS OF LAST FUNCTION PARAMETER........ MADR=IPAR(JFUNCT)
C     THE PARAMETERS OF THE FUNCTION ARE STORED IN IPAR(LADR) TO
C     IPAR(MADR).
C
C DATE: 1992, NOVEMBER 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:
      INTEGER NPAR
      PARAMETER (NPAR=8191)
      INTEGER IPAR(0:NPAR)
      REAL RPAR(0:NPAR)
      EQUIVALENCE (IPAR,RPAR)
      COMMON/VALC/IPAR
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: 1993, DECEMBER 18
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
      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     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
          READ(LUN,*) IPAR(LADR),IPAR(LADR+1),IPAR(LADR+2),RPAR(MADR)
          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))
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 STATAMENT - HERMITE REPRESENTATION:
*               CALL CURVN1(NX1,RPAR(JADR1),RPAR(JADR4),
*    *                      RPAR(JADR5),RPAR(MADR+1),SIGMA,IERR)
C             SECOND STATAMENT - 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)
      INTEGER LUN,NX1,NX2,NX3,IX1(NX1),IX2(NX2),IX3(NX3)
      REAL VAL(NX1,NX2,NX3)
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
C
      READ(LUN,*) (((VAL(IX1(I1),IX2(I2),IX3(I3)),I1=1,NX1),
     *                                            I2=1,NX2),I3=1,NX3)
      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 GEATER THAN OR
C             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:
      INTEGER NPAR
      PARAMETER (NPAR=8191)
      INTEGER IPAR(0:NPAR)
      REAL RPAR(0:NPAR)
      EQUIVALENCE (IPAR,RPAR)
      COMMON/VALC/IPAR
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: 1993, MAY 1
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
        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
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 STATAMENT - HERMITE REPRESENTATION:
*               CALL CURV2D(XX1,W(1),W(2),W(5),NX1,
*    *                      RPAR(JADR1),RPAR(JADR2),RPAR(JADR3),SIGMA)
C             SECOND STATAMENT - 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
