C SUBROUTINE FILE 'MODEL.FOR' TO SPECIFY A SEISMIC MODEL.
C
C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
C
C THIS FILE CONSISTS OF THE FOLLOWING EXTERNAL PROCEDURES:
C     MODELB..BLOCK DATA SUBROUTINE DEFINING COMMON BLOCKS MODELT AND
C             MODELC TO STORE THE DATA DESCRIBING THE MODEL.
C     MODEL1..SUBROUTINE DESIGNED TO READ THE INPUT DATA FOR THE
C             DESCRIPTION OF THE MODEL AND TO STORE THEM IN THE MEMORY.
C             SUBROUTINE MODEL1 READS THE INPUT DATA (1)-(8) LISTED
C             BELOW ITSELF, THEN CALLS SUBROUTINE SRFC1 (INCLUDED IN THE
C             SUBROUTINE FILE 'SRFC.FOR') TO READ THE INPUT DATA (9) FOR
C             SMOOTH SURFACES, AND FINALY CALLS SUBROUTINE PARM1
C             (INCLUDED IN THE SUBROUTINE FILE 'PARM.FOR') TO READ THE
C             INPUT DATA (10) FOR THE PARAMETERS OF THE MEDIUM.
C     NSRFC...INTEGER FUNCTION RETURNING THE NUMBER OF SURFACES COVERING
C             STRUCTURAL INTERFACES.
C     BLOCK...SUBROUTINE DESIGNED TO DETERMINE THE MUTUAL POSITION OF A
C             POINT AND A SIMPLE AND A COMPLEX BLOCK.
C     ISIDE...AUXILIARY FUNCTION TO SUBROUTINE BLOCK.
C     INTERF..AUXILIARY SUBROUTINE TO SUBROUTINE BLOCK.
C     VELOC...SUBROUTINE TRANSFORMING THE VALUES OF A MEDIUM PARAMETERS
C             INTO VELOCITIES AND LOSS FACTORS.
C     FPOWER..SUBROUTINE EVALUATING THE VALUE AND, POSSIBLY, THE THREE
C             FIRST AND SIX SECOND PARTIAL DERIVATIVES OF A FUNCTION, IF
C             THE VALUE AND THE THREE FIRST AND SIX SECOND PARTIAL
C             DERIVATIVES OF THE POWER-TH POWER OF THE FUNCTION ARE
C             ARE KNOWN.  PARTICULARLY, THIS IS AN AUXILIARY SUBROUTINE
C             TO THE SUBROUTINE VELOC.
C
C ATTENTION:
C     THE LINES DENOTED BY 'V' IN THE FIRST COLUMN SHOULD BE REMOVED
C     FROM THE SUBROUTINES VELOC AND FPOWER 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     VAR4 AND VAR5 OF THE FILE 'VAR.FOR' WOULD BE CALLED TO HANDLE THE
C     VARIATIONS.
C
C INPUT DATA FOR THE SPECIFICATION OF THE MODEL:
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 TEXTM), THE INPUT PARAMETER IS OF THE
C     TYPE REAL.
C (1) TEXTM
C     THE STRING CONTAINING THE NAME OF THE MODEL.  ONLY THE FIRST 80
C     CHARACTERS OF THE STRING ARE SIGNIFICANT.
C (2) KOORS,NEXPV,NEXPQ
C     KOORS...SPECIFIES THE TYPE OF THE RIGHT-HANDED COORDINATE SYSTEM
C             (SEE THE SUBROUTINE FILE 'METRIC.FOR'):
C             KOORS.LE.0: CARTESIAN COORDINATES.
C             KOORS.EQ.1: POLAR SPHERICAL COORDINATES (X1,X2,X3)=
C               (COLATITUDE,LONGITUDE,RADIUS).
C             KOORS.GE.2: GEOGRAPHIC SPHERICAL COORDINATES (X1,X2,X3)=
C               (LONGITUDE,LATITUDE,RADIUS).
C     NEXPV,NEXPQ... SPECIFY EXPONENTS OF THE POWER OF VELOCITIES
C             (NEXPV) AND LOSS FACTORS (NEXPQ) IN INPUT DATA.  FOR
C             EXAMPLE, UNIT VALUES OF NEXPV AND NEXPQ INDICATE THAT THE
C             PARAMETERS OF THE MEDIUM ARE VELOCITIES AND LOSS FACTORS,
C             INDICES EQUAL -1 INDICATE RECIPROCAL VALUES OF THESE
C             QUANTITIES, I.E.  SLOWNESSES AND QUALITY FACTORS.
C             WHILE USING THE BASIC SUBMITTED VERSION OF THE SUBROUTINE
C             FILE 'PARM.FOR', THE SPECIFICATION OF NEXPV=1, NEXPQ=1 IS
C             HIGHLY RECOMMENDED.
C     DEFAULT: KOORS=0, NEXPV=1, NEXPQ=1.
C (3) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX
C     BOUNDARIES OF THE MODEL.
C (4) NSRFC
C     NUMBER OF SMOOTH SURFACES IN THE MODEL.  THE SURFACES ARE INDEXED
C     SEQUENTIALLY IN ANY ORDER BY POSITIVE INTEGERS ISRFC FROM 1 TO
C     NSRFC.
C (5) NSB
C     NUMBER OF MATERIAL SIMPLE BLOCKS IN THE MODEL.  THE BLOCKS ARE
C     INDEXED SEQUENTIALLY IN ANY ORDER BY POSITIVE INTEGERS ISB FROM 1
C     TO NSB.  FREE-SPACE BLOCKS ARE NOT INDEXED.
C (6) NSB INPUT OPERATIONS (READ STATEMENTS):
C     FOR EACH SIMPLE BLOCK WITH INDEX ISB, THE INDICES OF THE SURFACES
C     FORMING THE SET F(+) AND THE INDICES OF THE SURFACES FORMING THE
C     SET F(-).  THE INDICES OF SURFACES FROM F(+) MUST BE POSITIVE, THE
C     INDICES OF SURFACES FROM F(-) MUST BE INDICATED BY NEGATIVE SIGNS.
C     THE INDICES MAY BE SPECIFIED IN AN ARBITRARY ORDER AND MUST BE
C     TERMINATED BY A SLASH.
C (7) NCB
C     NUMBER OF MATERIAL COMPLEX BLOCKS IN THE MODEL.  THE BLOCKS ARE
C     INDEXED SEQUENTIALLY BY POSITIVE INTEGERS ICB FROM 1 TO NCB.  THE
C     FREE-SPACE BLOCKS ARE NOT INDEXED.
C (8) NCB INPUT OPERATIONS (READ STATEMENTS):
C     FOR EACH COMPLEX BLOCK, THE INDICES OF SIMPLE BLOCKS FORMING THE
C     COMPLEX BLOCK.  THE INDICES MAY BE SPECIFIED IN AN ARBITRARY ORDER
C     AND MUST BE TERMINATED BY A SLASH.
C (9) THE DATA SPECIFYING NSRFC FUNCTIONS F(X1,X2,X3) DESCRIBING THE
C     SMOOTH SURFACES IN THE MODEL. THE DATA ARE READ BY SUBROUTINE
C     SRFC1.  FOR THEIR DESCRIPTION REFER TO SUBROUTINE SRFC1 (INCLUDED
C     IN THE SUBROUTINE FILE 'SRFC.FOR').
C (10) THE DATA SPECIFYING THE DISTRIBUTION OF PARAMETERS OF THE MODEL
C     IN ALL NCB MATERIAL COMPLEX BLOCKS.  THE DATA ARE READ BY
C     SUBROUTINE PARM1.  FOR THEIR DESCRIPTION REFER TO SUBROUTINE
C     PARM1 (INCLUDED IN THE SUBROUTINE FILE 'PARM.FOR').
C FOR AN EXAMPLE REFER TO THE SAMPLE INPUT DATA FOR THE MODEL.
C
C STORAGE IN THE MEMORY:
C     THE INPUT DATA (1) TO (8) DESCRIBING THE STRUCTURE OF THE MODEL
C     ARE STORED IN COMMON BLOCKS /MODELT/ AND /MODELC/ DEFINED IN THE
C     FOLLOWING SUBROUTINE:
C     ------------------------------------------------------------------
      BLOCK DATA MODELB
      CHARACTER*80 TEXTM
      COMMON/MODELT/TEXTM
      SAVE /MODELT/
      INTEGER MSB,MCB
      PARAMETER (MSB=128)
      PARAMETER (MCB=128)
      INTEGER NEXPV,NEXPQ
      REAL BOUNDM(6)
      INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB)
      EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB)
      COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB
      SAVE /MODELC/
      END
C     ------------------------------------------------------------------
C     TEXTM...THE NAME OF THE MODEL. STRING OF 80 CHARACTERS.
C     NEXPV,NEXPQ... SPECIFY EXPONENTS OF THE POWER OF VELOCITIES
C             (NEXPV) AND Q-FACTORS (NEXPQ) IN INPUT DATA.  FOR EXAMPLE,
C             UNIT VALUES OF NEXPV AND NEXPQ INDICATE THAT THE
C             PARAMETERS OF THE MEDIUM ARE VELOCITIES AND Q FACTORS,
C             INDICES EQUAL -1 INDICATE RECIPROCAL VALUES OF THESE
C             QUANTITIES, I.E. SLOWNESSES AND LOSS FACTORS.
C     BOUNDM..BOUNDARIES X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX OF THE
C             MODEL.
C     NSRFCS..NUMBER OF SMOOTH SURFACES IN THE MODEL. THE SURFACES
C             ARE INDEXED SEQUENTIALLY BY POSITIVE INTEGERS, FROM 1 TO
C             NSRFCS.  NSRFCS IS THE STORAGE LOCATION FOR NSRFC.
C     NSB...  NUMBER OF MATERIAL SIMPLE BLOCKS IN THE MODEL. THE BLOCKS
C             ARE INDEXED SEQUENTIALLY BY POSITIVE INTEGERS ISB FROM 1
C             TO NSB.  FREE-SPACE BLOCKS ARE NOT INDEXED.
C     KSB...  CONTAINS THE INDICES OF THE SURFACES BOUNDING INDIVIDUAL
C             SIMPLE BLOCKS.  KSB(ISB), FOR ISB = 1 TO NSB, SPECIFY THE
C             PARTITION OF ARRAY KSB(NSB+1:NSB+NS) AMONG THE SIMPLE
C             BLOCKS.  HERE NS IS THE TOTAL NUMBER OF ALL OCCURENCIES OF
C             THE INDICES OF THE SURFACES BOUNDING ALL INDIVIDUAL SIMPLE
C             BLOCKS IN THE INPUT DATA.  THE INDICES OF THE SURFACES
C             BOUNDING INDIVIDUAL SIMPLE BLOCKS ARE STORED FROM
C             KSB(NSB+1) TO KSB(NSB+NS). THE LOCATIONS KSB(NSB+NS+1:MSB)
C             ARE UNDEFINED.  IT MUST BE NSB+NS.LT.MSB.  THE INDICES OF
C             THE SURFACES BOUNDING THE SIMPLE BLOCK ISB ARE STORED IN
C             KSB(I1) TO KSB(I2), WITH
C               I1 = KSB(ISB-1)+1  ,
C               I2 = KSB(ISB)      ,
C             WHERE KSB(ISB-1)=NSB FOR ISB=1.  FOR EACH SIMPLE BLOCK
C             WITH INDEX ISB, THE INDICES OF THE SURFACES FORMING THE
C             SET F(+) ARE STORED WITH POSITIVE SIGNS, THE INDICES OF
C             SURFACES FROM F(-) WITH NEGATIVE SIGNS.  FOR AN EXAMPLE
C             REFER TO THE SAMPLE INPUT DATA FOR THE MODEL.
C     NCB...  NUMBER OF MATERIAL COMPLEX BLOCKS IN THE MODEL. THE BLOCKS
C             ARE INDEXED SEQUENTIALLY BY POSITIVE INTEGERS ICB FROM 1
C             TO NCB.  THE FREE-SPACE BLOCKS ARE NOT INDEXED.
C     KCB...  CONTAINS THE INDICES OF THE SIMPLE BLOCKS FORMING
C             INDIVIDUAL COMPLEX BLOCKS.  KCB(ICB), FOR ICB = 1 TO NCB,
C             SPECIFY THE PARTITION OF ARRAY KSB(NCB+1:NCB+NB) AMONG THE
C             COMPLEX BLOCKS.  HERE NB IS THE TOTAL NUMBER OF ALL
C             OCCURENCIES OF THE INDICES OF THE SIMPLE BLOCKS FORMING
C             INDIVIDUAL COMPLEX BLOCKS IN THE INPUT DATA.  THE INDICES
C             OF THE SIMPLE BLOCKS FORMING INDIVIDUAL COMPLEX BLOCKS ARE
C             STORED FROM KSB(NCB+1) TO KSB(NCB+NB). THE LOCATIONS
C             KSB(NCB+NB+1:MCB) ARE UNDEFINED. IT MUST BE NCB+NB.LT.MCB.
C             THE INDICES OF THE SIMPLE BLOCKS FORMING COMPLEX BLOCK ICB
C             ARE STORED IN KCB(I1) TO KCB(I2), WHERE
C               I1 = KCB(ICB-1)+1  ,
C               I2 = KCB(ICB)      .
C             HERE KCB(ICB-1)=NCB FOR ICB=1.  FOR AN EXAMPLE REFER TO
C             THE SAMPLE INPUT DATA FOR THE MODEL.
C     COMMON BLOCK /MODELT/ IS INCLUDED IN EXTERNAL PROCEDURE MODEL1.
C     COMMON BLOCK /MODELC/ IS INCLUDED IN EXTERNAL PROCEDURES MODEL1,
C     BLOCK, ISIDE, INTERF, AND MAY BE INCLUDED IN ANY OTHER SUBROUTINE.
C     ALL THE INPUT DATA ARE STORED SEQUENTIALLY IN THE SAME ORDER AS
C     THEY WERE READ.  THE ONLY EXCEPTION ARE LOCATIONS KSB(1) TO
C     KSB(NSB) AND KCB(1) TO KCB(NCB) WHICH ARE INSERTED WHEN THE INPUT
C     DATA ARE BEING READ .  THE INDEX OF THE LAST ALLOCATED NUMERIC
C     STORAGE UNIT OF ARRAY KSB IS NAMED MSB.  THE INDEX OF THE LAST
C     ALLOCATED NUMERIC STORAGE UNIT OF ARRAY KCB IS NAMED MCB.  THE
C     VALUES OF MSB AND MCB ARE GIVEN BY THE SIXTH AND SEVENTH STATEMENT
C     OF THE BLOCK DATA SUBROUTINE MODELB.  IF THE VALUE OF MSB OR MCB
C     IS CHANGED, IT MUST BE ADJUSTED IN ALL SUBROUTINES WHICH INCLUDE
C     THE COMMON BLOCK /MODELC/.
C
C DATE: 1992, NOVEMBER 30
C CODED BY LUDEK KLIMES
C
C=======================================================================
C
      SUBROUTINE MODEL1(LUN)
      INTEGER LUN
C
C SUBROUTINE MODEL1 READS THE INPUT DATA (1)-(8) FOR THE DESCRIPTION
C OF THE MODEL AND STORES THEM IN COMMON BLOCKS /MODELT/ AND /MODELC/.
C THEN IT CALLS SUBROUTINE SRFC1 (INCLUDED IN THE SUBROUTINE FILE
C 'SRFC.FOR') TO READ THE INPUT DATA (9) FOR SMOOTH SURFACES AND TO
C STORE THEM IN THE MEMORY.  FINALY, IT CALLS SUBROUTINE PARM1 (INCLUDED
C IN THE SUBROUTINE FILE 'PARM.FOR') TO READ THE INPUT DATA (10) FOR THE
C PARAMETERS OF THE MEDIUM AND TO STORE THEM IN THE MEMORY.
C
C INPUT:
C     LUN...  LOGICAL UNIT NUMBER OF THE EXTERNAL INPUT DEVICE
C             CONTAINING THE INPUT DATA.
C THE INPUT PARAMETER IS NOT ALTERED.
C
C NO OUTPUT.
C
C COMMON BLOCK /MODELT/:
      CHARACTER*80 TEXTM
      COMMON/MODELT/TEXTM
C COMMON BLOCK /MODELC/:
      INTEGER MSB,MCB
      PARAMETER (MSB=128)
      PARAMETER (MCB=128)
      INTEGER NEXPV,NEXPQ
      REAL BOUNDM(6)
      INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB)
      EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB)
      COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB
C ALL THE STORAGE LOCATIONS OF THE COMMON BLOCKS ARE DEFINED IN THIS
C SUBROUTINE.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL MODELB,METR1,SRFC1,PARM1
C     MODELB..BLOCK DATA SUBROUTINE OF THIS FILE.
C     METR1...FILE 'METRIC.FOR'.
C     SRFC1 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR').
C     PARM1 AND SUBSEQUENT ROUTINES... FILE 'PARM.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR').
C
C ERROR MESSAGES:
C     310...  INSUFFICIENT MEMORY IN /MODELC/:
C             INSUFFICIENT MEMORY FOR THE INPUT DATA IN COMMON BLOCK
C             /MODELC/.  THE DIMENSIONS MSB OR MCB OF ARRAYS KSB OR KCB,
C             RESPECTIVELY, MUST BE ENLARGED.  SEE THE BLOCK DATA
C             SUBROUTINE MODELB.
C     311...  BLOCK BOUNDED BY WRONG INTERFACE:
C             INDEX OF THE SURFACE BOUNDING THE SIMPLE BLOCK IS GREATER
C             THAN THE SPECIFIED NUMBER OF SURFACES.
C     312...  C. BLOCK COMPOSED OF WRONG S. BLOCK:
C             INDEX OF A SIMPLE BLOCK COMPOSING THE COMPLEX BLOCK IS
C             GREATER THAN THE SPECIFIED NUMBER OF SIMPLE BLOCKS.
C
C DATE: 1993, DECEMBER 18
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER I,J,L
C
      READ(LUN,*) TEXTM
      I=0
      NEXPV=1
      NEXPQ=1
      READ(LUN,*) I,NEXPV,NEXPQ
      CALL METR1(I)
      READ(LUN,*) BOUNDM
      READ(LUN,*) NSRFCS
C
C     SIMPLE BLOCKS:
C     NUMBER OF SIMPLE BLOCKS
      READ(LUN,*) NSB
C     INITIALIZING MEMORY FOR INDICES OF SURFACES BOUNDING SIMPLE BLOCKS
      L=NSB+1
      DO 11 I=L,MSB
        KSB(I)=0
   11 CONTINUE
C     READING INDICES OF SURFACES BOUNDING SIMPLE BLOCKS:
      DO 14 J=1,NSB
        READ(LUN,*) (KSB(I),I=L,MSB)
        DO 12 I=L,MSB
          IF(IABS(KSB(I)).GT.NSRFCS) THEN
            PAUSE 'ERROR 311 IN MODEL: BLOCK BOUNDED BY WRONG INTERFACE'
            STOP
          ELSE IF(KSB(I).EQ.0) THEN
            KSB(J)=I-1
            L=I
            GO TO 13
          END IF
   12   CONTINUE
        GO TO 99
   13   CONTINUE
   14 CONTINUE
C
C     COMPLEX BLOCKS:
C     NUMBER OF COMPLEX BLOCKS
      READ(LUN,*) NCB
C     INITIALIZING MEMORY FOR INDICES OF SIMPLE BLOCKS FORMING C. BLOCKS
      L=NCB+1
      DO 21 I=L,MCB
        KCB(I)=0
   21 CONTINUE
C     READING INDICES OF SIMPLE BLOCKS FORMING COMPLEX BLOCKS
      DO 24 J=1,NCB
        READ(LUN,*) (KCB(I),I=L,MCB)
        DO 22 I=L,MCB
          IF(KCB(I).LT.0.OR.KCB(I).GT.NSB) THEN
            STOP
     *        'ERROR 312 IN MODEL: C. BLOCK COMPOSED OF WRONG S. BLOCK.'
          ELSE IF(KCB(I).EQ.0) THEN
            KCB(J)=I-1
            L=I
            GO TO 23
          END IF
   22   CONTINUE
        GO TO 99
   23   CONTINUE
   24 CONTINUE
C
C     SMOOTH SURFACES:
      CALL SRFC1(LUN,NSRFCS)
C
C     MATERIAL PARAMETERS:
      CALL PARM1(LUN,NCB)
      RETURN
C
   99 CONTINUE
      PAUSE 'ERROR 310 IN MODEL1: INSUFFICIENT MEMORY IN /MODELC/'
      STOP
      END
C
C=======================================================================
C
      INTEGER FUNCTION NSRFC()
C
C INTEGER FUNCTION NSRFC IS DESIGNED TO RETURN THE NUMBER OF SURFACES
C COVERING STRUCTURAL INTERFACES.
C
C NO INPUT.
C
C OUTPUT:
C     NSRFC...NUMBER OF SURFACES COVERING STRUCTURAL INTERFACES.
C
C COMMON BLOCK /MODELC/:
      INTEGER MSB,MCB
      PARAMETER (MSB=128)
      PARAMETER (MCB=128)
      INTEGER NEXPV,NEXPQ
      REAL BOUNDM(6)
      INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB)
      EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB)
      COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED.
C
C DATE: 1989, DECEMBER 19
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     NO AUXILIARY STORAGE LOCATIONS.
C
      NSRFC=NSRFCS
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE BLOCK(COOR,ISRF1,ISB1,ISRF2,ISB2,ICB2,F)
      REAL COOR(3)
      INTEGER ISRF1,ISB1,ISB2,ICB2,ISRF2
      REAL F(10)
C
C THIS SUBROUTINE SEARCHES FOR THE SIMPLE BLOCK AND THE COMPLEX BLOCK IN
C WHICH A GIVEN POINT IS SITUATED.  THIS ROUTINE MAY BE ALSO USED TO
C DETERMINE THE INDEX OF A BLOCK TOUCHING A SPECIFIED BLOCK AT A GIVEN
C POINT SITUATED ON THE BOUNDARY OF THE SPECIFIED BLOCK (THE SITUATION
C WHICH MAY OCCUR WHEN A RAY IMPINGES ON A BOUNDARY OF A BLOCK).
C ANOTHER FUNCTION OF THE ROUTINE IS TO DETERMINE THE INDEX OF ANY OF
C THE SURFACES BOUNDING A BLOCK AND SEPARATING THE BLOCK FROM THE GIVEN
C POINT.
C
C INPUT:
C     COOR... ARRAY CONTAINING COORDINATES X1, X2, X3 OF A GIVEN POINT.
C     ISRF1...ISRF1.NE.0: INDEX OF THE SURFACE AT WHICH THE GIVEN POINT
C               IS SITUATED.  THE SIGN IS IGNORED.
C             ISRF1.EQ.0: THE GIVEN POINT IS NOT SITUATED AT ANY SURFACE
C     ISB1... FOR ISRF1.NE.0: INDEX OF A SIMPLE BLOCK BOUNDED BY THE
C               SURFACE IABS(ISRF1) - SEARCH FOR THE INDEX OF A
C               NEIGHBOURING SIMPLE BLOCK TOUCHING ISB1 AT THE GIVEN
C               POINT.
C             FOR ISRF1.EQ.0: INDEX OF AN ARBITRARY SIMPLE BLOCK OR
C               ISB1=0.
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C OUTPUT:
C     ISRF2...FOR THE GIVEN POINT NOT SITUATED INSIDE THE BLOCK ISB1 OR
C             AT ITS BOUNDARY, ISRF2 HAS THE MEANING OF THE INDEX OF ONE
C             OF THE SURFACES BOUNDING THE SIMPLE BLOCK ISB1 AND
C             SEPARATING THE GIVEN POINT FROM THE SIMPLE BLOCK ISB1,
C             SUPPLEMENTED BY A SIGN '+' OR '-' FOR THE SIMPLE BLOCK
C             ISB1 SITUATED AT THE POSITIVE OR NEGATIVE SIDE OF THE
C             SURFACE, RESPECTIVELY.
C             OTHERWISE ZERO.
C     ISB2... FOR ISRF1.NE.0, ISRF2.EQ.0: INDEX OF A SIMPLE BLOCK
C               NEIGHBOURING TO THE SIMPLE BLOCK ISB1 AND SEPARATED FROM
C               THE SIMPLE BLOCK ISB1 BY SURFACE IABS(ISRF1) AT THE
C               GIVEN POINT.
C             FOR ISRF1.NE.0, ISRF2.NE.0: INDEX OF A SIMPLE BLOCK IN
C               WHICH THE GIVEN POINT IS SITUATED.  IN THIS CASE, THE
C               GIVEN POINT MAY BE SITUATED EITHER INSIDE THE SIMPLE
C               BLOCK ISB2 OR IN THE VICINITY OF ITS BOUNDARY FORMED BY
C               THE SURFACE IABS(ISRF1).  FROM THE TWO POSSIBLE SIMPLE
C               BLOCKS TOUCHING THE SURFACE IABS(ISRF1) AT THE GIVEN
C               POINT, THAT BEING SITUATED AT THE SAME SIDE OF THE
C               SURFACE ISRF1 AS THE SIMPLE BLOCK ISB1, IS SELECTED.
C             FOR ISRF1.EQ.0: INDEX OF THE SIMPLE BLOCK IN WHICH THE
C               GIVEN POINT IS SITUATED.
C             IF THERE IS NO SIMPLE MATERIAL BLOCK OF THE ABOVE
C             PROPERTIES, ISB2=0.
C     ICB2... INDEX OF THE COMPLEX BLOCK IN WHICH THE SIMPLE BLOCK ISB2
C             IS SITUATED.  ICB2=0 IF THE SIMPLE BLOCK ISB2 IS NOT
C             SITUATED IN ANY COMPLEX BLOCK.  FOR ISB2=0: ICB2=0.
C     F...    FOR ISRF2.NE.0: ARRAY CONTAINING THE VALUE AND PARTIAL
C               DERIVATIVES F, F1, F3, F11, F12, F22, F13, F23, F33 OF
C               THE FUNCTION DESCRIBING THE SURFACE IABS(ISRF2).
C             FOR ISRF2.EQ.0: UNDEFINED.
C
C COMMON BLOCK /MODELC/:
      INTEGER MSB,MCB
      PARAMETER (MSB=128)
      PARAMETER (MCB=128)
      INTEGER NEXPV,NEXPQ
      REAL BOUNDM(6)
      INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB)
      EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB)
      COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL ISIDE,INTERF,SRFC2
      INTEGER ISIDE
C     ISIDE,INTERF... THIS FILE.
C     SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT
C              FILES (LIKE 'VAL.FOR' AND 'FIT.FOR').
C
C ERROR MESSAGES:
C     313...  WRONG INDEX OF SURFACE:
C             ABSOLUTE VALUE OF THE INPUT PARAMETER ISRFC1 (INDEX OF THE
C             SURFACE) IS GREATER THAN THE NUMBER NSRFC OF THE SURFACES
C             COVERING STRUCTURAL INTERFACES.
C     314...  WRONG INDEX OF SIMPLE BLOCK:
C             PARAMETER ISB1 (INDEX OF THE SIMPLE BLOCK) IS EITHER
C             NEGATIVE OR GREATER THAN THE NUMBER NSB OF SIMPLE BLOCKS.
C
C DATE: 1993, MAY 18
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER I,J,I1,I2,ISIDE1,ISIDE2,ISRF
      REAL FAUX(10)
C
C     CHECKING INPUT VALUES:
      IF(ISRF1.LT.-NSRFCS.OR.ISRF1.GT.NSRFCS) THEN
        PAUSE 'ERROR 313 IN BLOCK: WRONG INDEX OF SURFACE'
        STOP
      END IF
      IF(ISB1.LT.0.OR.ISB1.GT.NSB) THEN
        PAUSE 'ERROR 314 IN BLOCK: WRONG INDEX OF SIMPLE BLOCK'
        STOP
      END IF
C
C     POSITION OF THE GIVEN POINT WITH RESPECT TO THE GIVEN SIMPLE BLOCK
C     ISB1:
      IF(ISB1.EQ.0) THEN
C       NO SIMPLE BLOCK SPECIFIED:
        ISRF2=0
      ELSE
        CALL INTERF(COOR,ISRF1,ISB1,ISRF2,F)
        IF(ISRF1.EQ.0) THEN
          IF(ISRF2.EQ.0) THEN
C           THE POINT IS INSIDE SIMPLE BLOCK ISB1:
            ISB2=ISB1
            GO TO 10
          END IF
        ELSE
          IF(ISRF2.NE.0) THEN
C           THE POINT BEING SITUATED AT SURFACE ISRF1 BOUNDING SIMPLE
C           BLOCK ISB1 IS NOT SITUATED AT THE BOUNDARY OF SIMPLE BLOCK
C           ISB1:
            ISIDE1=ISIDE(ISRF1,ISB1)
          ELSE
C           THE POINT IS SITUATED AT THE BOUNDARY OF SIMPLE BLOCK ISB1:
            ISIDE1=-ISIDE(ISRF1,ISB1)
          END IF
        END IF
      END IF
C
C     SEARCH FOR THE SIMPLE BLOCK IN WHICH THE GIVEN POINT IS
C     SITUATED:
      I2=MAX0(ISB1-1,NSB-ISB1)
      DO 2 J=1,I2
        DO 1 I=1,-1,-2
          ISB2=ISB1+I*J
          IF(ISB2.GT.0.AND.ISB2.LE.NSB) THEN
            IF(ISRF1.NE.0) THEN
              ISIDE2=ISIDE(ISRF1,ISB2)
            END IF
            IF(ISRF1.EQ.0.OR.ISIDE1.EQ.ISIDE2) THEN
              CALL INTERF(COOR,ISRF1,ISB2,ISRF,FAUX)
              IF(ISRF.EQ.0) GO TO 10
            END IF
          END IF
    1   CONTINUE
    2 CONTINUE
C     NO SIMPLE BLOCK HAS BEEN FOUND:
      ISB2=0
      ICB2=0
      RETURN
C
C     DETERMINATION OF THE COMPLEX BLOCK:
   10 CONTINUE
      DO 12 J=1,NCB
        I1=KCB(J-1)+1
        I2=KCB(J)
        DO 11 I=I1,I2
          IF(KCB(I).EQ.ISB2) THEN
            ICB2=J
            RETURN
          END IF
   11   CONTINUE
   12 CONTINUE
C     NO COMPLEX BLOCK:
      ICB2=0
C
      RETURN
      END
C
C=======================================================================
C
      INTEGER FUNCTION ISIDE(ISRF,ISB)
      INTEGER ISRF,ISB
C
C THIS IS AN AUXILIARY FUNCTION TO THE SUBROUTINE BLOCK.
C THIS FUNCTION DETERMINES THE MUTUAL POSITION OF A SURFACE AND A SIMPLE
C BLOCK.
C
C INPUT:
C     ISRF... INDEX OF THE SURFACE.  THE SIGN IS IGNORED.
C     ISB...  INDEX OF THE SIMPLE BLOCK.
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C OUTPUT:
C     ISIDE...ISIDE=-1: THE SIMPLE BLOCK IS BOUNDED BY THE SURFACE AND
C               IS SITUATED ON ITS NEGATIVE SIDE.
C             ISIDE= 1: THE SIMPLE BLOCK IS BOUNDED BY THE SURFACE AND
C               IS SITUATED ON ITS POSITIVE SIDE.
C             ISIDE= 2: THE SIMPLE BLOCK IS NOT BOUNDED BY THE SURFACE.
C
C COMMON BLOCK /MODELC/:
      INTEGER MSB,MCB
      PARAMETER (MSB=128)
      PARAMETER (MCB=128)
      INTEGER NEXPV,NEXPQ
      REAL BOUNDM(6)
      INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB)
      EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB)
      COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED.
C
C DATE: 1989, DECEMBER 19
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER IS,LS,MS
C
      LS=KSB(ISB-1)+1
      MS=KSB(ISB)
C
C     LOOP FOR SURFACES BOUNDING SIMPLE BLOCK ISB:
      DO 1 IS=LS,MS
        IF(IABS(KSB(IS)).EQ.IABS(ISRF)) THEN
          ISIDE=ISIGN(1,KSB(IS))
          RETURN
        END IF
    1 CONTINUE
C
      ISIDE=2
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE INTERF(COOR,ISRF1,ISB,ISRF2,F)
      REAL COOR(3),F(10)
      INTEGER ISRF1,ISB,ISRF2
C
C THIS IS AN AUXILIARY SUBROUTINE TO THE SUBROUTINE BLOCK.
C THIS SUBROUTINE DETERMINES THE POSITION OF A GIVEN POINT WITH RESPECT
C TO A GIVEN SIMPLE BLOCK.
C
C INPUT:
C     COOR... ARRAY CONTAINING COORDINATES X1, X2, X3 OF A GIVEN POINT.
C     ISRF1...ISRF1.NE.0: INDEX OF THE SURFACE AT WHICH THE GIVEN POINT
C               IS SITUATED.  THE SIGN IS IGNORED.
C             ISRF1.EQ.0: THE GIVEN POINT IS NOT SITUATED AT ANY
C               SURFACE.
C     ISB...  INDEX OF THE GIVEN SIMPLE BLOCK.
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C OUTPUT:
C     ISRF2...INDEX OF A SURFACE SEPARATING THE GIVEN POINT AND THE
C               SIMPLE BLOCK ISB, SUPPLEMENTED BY A SIGN '+' OR '-' FOR
C               THE SIMPLE BLOCK ISB1 SITUATED AT THE POSITIVE OR
C               NEGATIVE SIDE OF THE SURFACE, RESPECTIVELY.
C             ISRF2=0  IF THE GIVEN POINT IS SITUATED INSIDE THE SIMPLE
C               BLOCK ISB.
C     F...    FOR ISRF2.NE.0: ARRAY CONTAINING THE VALUE AND PARTIAL
C               DERIVATIVES F, F1, F3, F11, F12, F22, F13, F23, F33 OF
C               THE FUNCTION DESCRIBING THE SURFACE IABS(ISRF2) AT THE
C               GIVEN POINT.
C             FOR ISRF2.EQ.0: UNDEFINED.
C
C COMMON BLOCK /MODELC/:
      INTEGER MSB,MCB
      PARAMETER (MSB=128)
      PARAMETER (MCB=128)
      INTEGER NEXPV,NEXPQ
      REAL BOUNDM(6)
      INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB)
      EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB)
      COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL SRFC2
C     SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT
C              FILES (LIKE 'VAL.FOR' AND 'FIT.FOR').
C
C DATE: 1989, DECEMBER 19
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER IS,JS,KS,LS,MS
C
      LS=KSB(ISB-1)+1
      MS=KSB(ISB)
C
C     LOOP FOR SURFACES BOUNDING SIMPLE BLOCK ISB:
      DO 1 IS=LS,MS
        KS=KSB(IS)
        JS=IABS(KS)
        IF(JS.NE.IABS(ISRF1)) THEN
          CALL SRFC2(JS,COOR,F)
          IF(F(1)*FLOAT(KS).LT.0.) THEN
            ISRF2=KS
            RETURN
          END IF
        END IF
    1 CONTINUE
C
      ISRF2=0
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE VELOC(IWAVE,UP,US,QP,QS,VP,VS,VD,QL)
      INTEGER IWAVE
      REAL UP(10),US(10),QP,QS,VP,VS,VD(10),QL
C
C THIS SUBROUTINE TRANSFORMS THE VALUES OF PARAMETERS OF THE MEDIUM INTO
C VELOCITIES AND LOSS FACTORS.
C
C INPUT:
C     IWAVE...TYPE OF WAVE.
C             IWAVE.GE.0: P WAVE,
C             IWAVE.LT.0: S WAVE.
C     UP,US...POWERS OF P AND S WAVE VELOCITIES AND THEIR FIRST AND
C             SECOND PARTIAL DERIVATIVES (THE EXPONENT OF THE POWERS IS
C             NEXPV, SEE 'INPUT DATA FOR THE MODEL'), IN ORDER U, U1,
C             U2, U3, U11, U12, U22, U13, U23, U33.
C     QP,QS...POWERS OF THE LOSS FACTORS OF P AND S WAVES (THE EXPONENT
C             OF THE POWERS IS NEXPQ, SEE 'INPUT DATA FOR THE MODEL').
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C OUTPUT:
C     VP,VS...P AND S WAVE VELOCITIES.
C     VD...   VELOCITY AND ITS FIRST AND SECOND PARTIAL DERIVATIVES
C             ORDERED AS UP, US, CORRESPONDING TO THE WAVE SPECIFIED BY
C             IWAVE, IN ORDER V, V1, V2, V3, V11, V12, V22, V13, V23,
C             V33.
C     QL...   LOSS FACTOR CORRESPONDING TO THE WAVE SPECIFIED BY IWAVE.
C
C COMMON BLOCK /MODELC/:
      INTEGER MSB,MCB
      PARAMETER (MSB=128)
      PARAMETER (MCB=128)
      INTEGER NEXPV,NEXPQ
      REAL BOUNDM(6)
      INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB)
      EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB)
      COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL FPOWER
C     FPOWER...THIS FILE.
C
C DATE: 1992, DECEMBER 31
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATION:
      REAL POWER,AUX1(1),AUX2(1)
C
      POWER=FLOAT(NEXPV)
      IF(IWAVE.GE.0) THEN
        CALL FPOWER(10,UP,POWER,VD)
V       CALL VAR5(1,1)
        VP=VD(1)
        CALL FPOWER(1,US,POWER,AUX2)
        VS=AUX2(1)
        AUX1(1)=QP
      ELSE
        CALL FPOWER(1,UP,POWER,AUX2)
        VP=AUX2(1)
        CALL FPOWER(10,US,POWER,VD)
V       CALL VAR5(2,2)
        VS=VD(1)
        AUX1(1)=QS
      END IF
      CALL FPOWER(1,AUX1,FLOAT(NEXPQ),AUX2)
      QL=AUX2(1)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE FPOWER(N,FINP,POWER,FOUT)
      INTEGER N
      REAL FINP(N),POWER,FOUT(N)
C
C THIS SUBROUTINE EVALUATES THE VALUE AND, POSSIBLY, THE THREE FIRST AND
C SIX SECOND PARTIAL DERIVATIVES OF A FUNCTION IF THE VALUE AND THE
C THREE FIRST AND SIX SECOND PARTIAL DERIVATIVES OF THE POWER-TH POWER
C OF THE FUNCTION ARE KNOWN.
C
C INPUT:
C     N...    FOR N=1: ONLY THE FUNCTION VALUE IS EVALUATED. THE
C               DERIVATIVES ARE IGNORED.
C             FOR N=4: THE VALUE AND THE THREE FIRST PARTIAL DERIVATIVES
C               ARE EVALUATED.
C             FOR N=10: THE VALUE AND THE THREE FIRST AND SIX SECOND
C               PARTIAL DERIVATIVES ARE EVALUATED.
C     FINP... ARRAY CONTAINING THE VALUE, THE FIRST AND SECOND PARTIAL
C             DERIVATIVES OF THE POWER-TH POWER OF THE FUNCTION TO BE
C             EVALUATED, IN THE ORDER F, F1, F2, F3, F11, F12, F22, F13,
C             F23, F33.  FOR N=1, ONLY THE FUNCTION VALUE IS REQUIRED.
C     POWER...THE SPECIFIED FUNCTION IS EQUAL TO THE POWER-TH POWER OF
C             THE CORRESPONDING PHYSICAL QUANTITY.
C NONE OF THE INPUT PARAMETERS ARE ALTERED (EXCEPT FINP IF THIS
C PARAMETER AND FOUT ARE IDENTICAL IN THE CALLING SEQUENCE).
C
C OUTPUT:
C     FOUT... ARRAY CONTAINING THE VALUE, THE FIRST AND SECOND PARTIAL
C             DERIVATIVES OF THE EVALUATED FUNCTION, IN THE ORDER F, F1,
C             F2, F3, F11, F12, F22, F13, F23, F33.  THIS PARAMETER MAY
C             COINCIDE WITH FINP, IN WHICH CASE FINP IS DESTROYED ON
C             OUTPUT.  NOTE THAT THIS COINCIDENCE IS AN EXCEPTION TO
C             ANSI STANDARD OF FORTRAN 77.
C
C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED.
C
C DATE: 1992, NOVEMBER 30
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER I
      REAL F,AUX1,AUX2
C
      IF(0.999.LT.POWER.AND.POWER.LT.1.001) THEN
        DO 1 I=1,N
          FOUT(I)=FINP(I)
    1   CONTINUE
V       CALL VAR4(0,1.)
      ELSE
        IF(-1.001.LT.POWER.AND.POWER.LT.-0.999) THEN
          F=1./FINP(1)
        ELSE
          F=FINP(1)**(1./POWER)
        END IF
        FOUT(1)=F
        IF(N.GT.1) THEN
          AUX1= F/(FINP(1)*POWER)
          AUX2= (POWER-1.)/F
          FOUT(2)=AUX1*FINP(2)
          FOUT(3)=AUX1*FINP(3)
          FOUT(4)=AUX1*FINP(4)
          IF(N.GT.4) THEN
            FOUT(5)=AUX1*FINP(5)-AUX2*FOUT(2)*FOUT(2)
            FOUT(6)=AUX1*FINP(6)-AUX2*FOUT(2)*FOUT(3)
            FOUT(7)=AUX1*FINP(7)-AUX2*FOUT(3)*FOUT(3)
            FOUT(8)=AUX1*FINP(8)-AUX2*FOUT(2)*FOUT(4)
            FOUT(9)=AUX1*FINP(9)-AUX2*FOUT(3)*FOUT(4)
            FOUT(10)=AUX1*FINP(10)-AUX2*FOUT(4)*FOUT(4)
          END IF
V         CALL VAR4(0,AUX1)
V         CALL VAR4(2,-AUX2*FOUT(2))
V         CALL VAR4(3,-AUX2*FOUT(3))
V         CALL VAR4(4,-AUX2*FOUT(4))
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
