C SUBROUTINE FILE 'CODE.FOR' - CODES FOR ELEMENTARY WAVES.
C
C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
C
C THIS FILE CONSISTS OF:
C     GENERAL DESCRIPTION OF THE CODES FOR ELEMENTARY WAVES.
C     CODEB...BLOCK DATA SUBROUTINE DEFINING COMMON BLOCK /COD/ TO STORE
C             THE CODE OF THE COMPUTED ELEMENTARY WAVE.
C     CODE1...SUBROUTINE DESIGNED TO READ THE INPUT DATA FOR THE CODES
C             OF ELEMENTARY WAVES AND TO STORE THEM IN THE COMMON BLOCK
C             /COD/.
C     LCODE...INTEGER FUNCTION RETURNING THE LENGTH OF THE CODE OF THE
C             CURRENT ELEMENTARY WAVE.
C     NELEM...INTEGER FUNCTION RETURNING THE NUMBER OF RAY ELEMENTS
C             CORRESPONDING TO THE GIVEN POSITION IN THE CODE OF THE
C             ELEMENTARY WAVE.
C     CODE... SUBROUTINE DESIGNED TO TRANSFORM THE USED NUMERICAL CODE
C             OF THE ELEMENTARY WAVE UNDER CONSIDERATION INTO
C             INSTRUCTIONS SPECIFYING THE BEHAVIOUR OF THE WAVE AT THE
C             INITIAL POINTS OF RAYS AND AT ALL POINTS OF INCIDENCE OF
C             THE RAYS AT INTERFACES (BOUNDARIES OF COMPLEX BLOCKS).
C
C INPUT DATA FOR THE CODES OF ELEMENTARY WAVES:
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 TEXTC), THE INPUT PARAMETER IS OF THE
C     TYPE REAL.
C (1) TEXTC
C     STRING DESCRIBING THE DATA.  ONLY THE FIRST 80 CHARACTERS OF THE
C     STRING ARE SIGNIFICANT.
C (2) NSKIP
C     THE NUMBER OF ELEMENTARY WAVES TO BE SKIPPED OVER.  THE
C     (NSKIP+1)-TH ELEMENTARY WAVE IN THESE INPUT DATA WILL BE THE FIRST
C     COMPUTED WAVE.  DEFAULT VALUE: 0.
C FOR EACH ELEMENTARY WAVE IWAVE THE FOLLOWING DATA (3):
C (3) KODTYP,KODE
C     KODTYP..DETERMINES THE TYPE OF THE CODE
C             1... BLOCK-SURFACE CODE,
C             2... BASIC BLOCK CODE,
C             3... COMPOUND-ELEMENT BLOCK CODE,
C             4... GENERALIZED LAYER CODE,
C             5... TRANSMISSION-REFLECTION CODE.
C     KODE... THE CODE OF THE ELEMENTARY WAVE IN THE FORM OF A FINITE
C             SEQUENCE OF INTEGERS TERMINATED BY A SLASH.
C     ATTENTION:
C             GENERALLY, DIFFERENT TYPES OF THE CODE MAY BE COMBINED
C             WITHIN THE SAME INPUT DATA, BUT THE USER HAS THEN BE VERY
C             CAREFUL.  IF COMBINING THE DIFFERENT TYPES OF THE CODE,
C             OUTPUT PARAMETERS IWAVE0 AND IKODE OF SUBROUTINE CODE1 MAY
C             HAVE NO SENSE.  THEN:
C             IF THE PARAMETER MODCRT OF THE DATA SET 'DCRT' (SEE SOURCE
C               CODE FILE 'RAY.FOR') IS ZERO, THE OUTPUT PARAMETERS
C               IWAVE0 AND IKODE OF SUBROUTINE CODE1 ARE HOPEFULLY
C               UNUSED.
C             IF THE PARAMETER MODCRT OF THE DATA SET 'DCRT' IS NOT
C               ZERO, MIXING THE TYPES OF THE WAVE CODE SHOULD BE
C               AVOIDED UNLESS THE USER FINDS A WAY HOW TO DO IT
C               REASONABLY.
C (4) A SLASH.
C
C STORAGE IN THE MEMORY:
C     THE INPUT DATA (2) FOR THE COMPUTED ELEMENTARY WAVE ARE STORED IN
C     THE COMMON BLOCK /COD/ DEFINED IN THE FOLLOWING SUBROUTINE:
C     ------------------------------------------------------------------
      BLOCK DATA CODEB
      INTEGER MKODE
      PARAMETER (MKODE=128)
      INTEGER KODTYP,KODE(MKODE)
      COMMON /COD/ KODTYP,KODE
      SAVE /COD/
      END
C     ------------------------------------------------------------------
C     KODTYP..DETERMINES THE TYPE OF THE CODE, SEE THE INPUT DATA (2).
C     KODE... ARRAY CONTAINING THE CODE OF THE ELEMENTARY WAVE.
C             THE CODE IS A FINITE SEQUENCE OF NONZERO INTEGERS.
C             ZERO INDICATES THE END OF A CODE.
C     COMMON BLOCK /COD/ IS INCLUDED IN EXTERNAL PROCEDURES COD1, CODE,
C     WRIT1,  AND MAY BE INCLUDED IN ANY OTHER SUBROUTINE. IF MKODE IS
C     CHANGED, IT MUST BE ADJUSTED IN ALL SUBROUTINES WHICH INCLUDE
C     COMMON BLOCK /COD/.
C
C DATE: 1993, AUGUST 3
C CODED BY LUDEK KLIMES
C
C=======================================================================
C
C CODES FOR ELEMENTARY WAVES - GENERAL DESCRIPTION
C
C     WE CONSIDER ORDINARY SEISMIC BODY WAVES, SUCH AS REFRACTED,
C     PRIMARILY OR MULTIPLY REFLECTED, POSSIBLY CONVERTED WAVES.  IN
C     GENERAL, INCIDENCE OF A WAVE AT AN INTERFACE (BOUNDARY OF A
C     COMPLEX BLOCK) PRODUCES FOUR WAVES, REFLECTED P AND S, AND
C     TRANSMITTED P AND S WAVES.  WHEN PERFORMING COMPLETE RAY TRACING,
C     WE MUST KNOW A PRIORI WHICH OF THE FOUR GENERATED WAVES TO FOLLOW.
C     THIS DECISION MUST BE MADE AT EACH INTERFACE.  THE ALPHANUMERIC
C     STRING SPECIFYING THE BEHAVIOUR OF A RAY FROM ITS INITIAL POINT TO
C     ITS ENDPOINT IS A 'CODE'.  IN THE FOLLOWING, WE SHALL CONSIDER THE
C     CODE TO BE A SEQUENCE OF NONZERO INTEGERS.
C
C     THE TERM 'ELEMENTARY WAVE' DOES NOT HAVE UNIQUE MEANING IN THE
C     LITERATURE.  HERE WE APPLY THE TERM ELEMENTARY WAVE TO THAT PART
C     OF THE WAVEFIELD THAT IS DESCRIBED BY ONE SPECIFIC CODE.  SINCE
C     THERE MAY BE VARIOUS TYPES OF CODES, THERE IS ALSO A VARIETY OF
C     DIVISIONS OF THE WAVE FIELD INTO ELEMENTARY WAVES.
C
C     WE INTRODUCE THE TERM 'ELEMENT OF A RAY', WHICH HAS AN IMPORTANT
C     MEANING IN THE CONSTRUCTION OF CODES.  BY AN ELEMENT OF A RAY, WE
C     DENOTE THAT PART OF THE RAY THAT IS SITUATED IN ONE COMPLEX BLOCK
C     BETWEEN TWO SUCCESSIVE POINTS OF REFLECTION/TRANSMISSION, OR
C     BETWEEN THE INITIAL POINT OR ENDPOINT OF THE RAY AND THE CLOSEST
C     POINT OF REFLECTION/TRANSMISSION, OR BETWEEN THE INITIAL POINT AND
C     THE ENDPOINT OF THE RAY, IF THE RAY IS ENTIRELY SITUATED IN ONE
C     COMPLEX BLOCK.
C
C     IN THE FOLLOWING, WE PRESENT SEVERAL POSSIBLE TYPES OF NUMERICAL
C     CODES OF ELEMENTARY WAVES.
C
C EXAMPLES OF CODE TYPES
C     (1) BLOCK-SURFACE CODE
C             TWO INTEGERS ARE USED FOR ANY ELEMENT OF THE RAY.  THE
C             ABSOLUTE VALUE OF THE FIRST INTEGER SPECIFIES THE INDEX OF
C             THE COMPLEX BLOCK, IN WHICH THE ELEMENT OF THE RAY IS
C             SITUATED, THE SIGN SPECIFIES THE TYPE OF THE WAVE (PLUS
C             SIGN...  P WAVE, MINUS SIGN...  S WAVE).  THE SECOND
C             NUMBER IS THE INDEX OF THE SMOOTH SURFACE, AT WHICH THE
C             ELEMENT TERMINATES.  THE CODE OF THE RAY IS A CHAIN OF THE
C             ABOVE DOUBLETS DESCRIBING THE RAY FROM ITS INITIAL POINT
C             TO ITS ENDPOINT.  THE INDEX OF A SURFACE, AT WHICH THE
C             LAST ELEMENT OF A RAY TERMINATES, NEED NOT BE SPECIFIED.
C             IN THIS CASE THE RAY IS ALLOWED TO TERMINATE ON ANY
C             SURFACE BOUNDING THE COMPLEX BLOCK IN WHICH THE LAST
C             ELEMENT OF THE RAY IS SITUATED.  APPLICATION OF THIS CODE
C             IS STRAIGHTFORWARD.
C
C             IF ANY COMPLEX BLOCK IS BOUNDED BY SEVERAL SMOOTH
C             SURFACES, THE BLOCK-SURFACE CODE MAY DIVIDE ARTIFICIALLY
C             THE WAVE FIELD INTO SEVERAL ELEMENTARY WAVES.  THIS WOULD
C             RESULT IN REPEATED APPLICATION OF THE COMPLETE RAY TRACING
C             TO ALL OF THE ELEMENTARY WAVES, WHAT MAY BE TIME
C             CONSUMING.  IN SUCH SITUATIONS, IT IS THEREFORE DESIRABLE
C             TO USE OTHER CODES, FOR WHICH THE ELEMENTARY WAVE HAS A
C             MORE GENERAL MEANING.
C
C     (2) BASIC BLOCK CODE
C             THIS CODE MAY BE OBTAINED FROM THE BLOCK-SURFACE CODE IF
C             THE INDICES OF THE SURFACES, AT WHICH INDIVIDUAL ELEMENTS
C             OF A RAY TERMINATE, ARE OMITTED.  THEN, ANY ELEMENT OF THE
C             RAY IS DESCRIBED BY A SINGLE INTEGER SPECIFYING THE
C             COMPLEX BLOCK AND WAVE TYPE, AND IT MAY TERMINATE AT ANY
C             SURFACE BOUNDING THE COMPLEX BLOCK.  THE CODE OF THE RAY
C             IS AGAIN A CHAIN OF THE NUMBERS DESCRIBING SUCCESIVELY ITS
C             INDIVIDUAL ELEMENTS.
C
C             ELEMENTARY WAVES SPECIFIED BY THE BLOCK-SURFACE CODE,
C             PASSING THROUGH THE SAME COMPLEX BLOCKS, AND REFLECTED
C             FROM DIFFERENT SURFACES BOUNDING A COMPLEX BLOCK, ARE
C             UNITED INTO ONE ELEMENTARY WAVE IN THE BASIC BLOCK CODE.
C             SIMILARLY, ELEMENTARY WAVES SPECIFIED BY THE BLOCK-SURFACE
C             CODE AND TRANSMITTED FROM A BLOCK INTO THE NEIGHBOURING
C             BLOCK ACROSS DIFFERENT SURFACES SEPARATING THESE BLOCKS,
C             ARE UNITED INTO ONE ELEMENTARY WAVE BY THE BASIC BLOCK
C             CODE.
C
C     (3) COMPOUND-ELEMENT BLOCK CODE
C             LET US INTRODUCE THE TERMS SIMPLE AND COMPOUND ELEMENT OF
C             A RAY.  FOR THIS PURPOSE, WE CALL THE 'LOWER-INDEX
C             BOUNDARY' OF A COMPLEX BLOCK THAT PART OF ITS BOUNDARY,
C             WHICH SEPARATES THE COMPLEX BLOCK EITHER FROM COMPLEX
C             BLOCKS WITH LOWER INDICES OR FROM A FREE SPACE.  THE
C             REMAINING PART OF THE BOUNDARY IS CALLED THE
C             'HIGHER-INDEX BOUNDARY' OF THE COMPLEX BLOCK.  THE INITIAL
C             POINT OF A RAY IS TREATED IN THE SAME WAY AS THE POINTS
C             SITUATED AT A LOWER-INDEX BOUNDARY.
C
C             AN ELEMENT OF A RAY IS CALLED 'SIMPLE ELEMENT' IF ITS
C             INITIAL POINT (E.G. THE POINT WHERE THE WAVE ENTERS THE
C             CORRESPONDING COMPLEX BLOCK) AND ITS ENDPOINT (E.G. THE
C             POINT AT WHICH THE WAVE LEAVES THE COMPLEX BLOCK) ARE
C             SITUATED ONE ON THE LOWER-INDEX AND THE OTHER ON THE
C             HIGHER-INDEX BOUNDARY OF THE COMPLEX BLOCK.  THE
C             'COMPOUND ELEMENT' IS SUCH AN ELEMENT FOR WHICH ITS
C             INITIAL POINT AND ITS ENDPOINT ARE BOTH SITUATED EITHER
C             ON THE LOWER-INDEX OR ON THE HIGHER-INDEX BOUNDARY OF THE
C             COMPLEX BLOCK.  THE COMPOUND ELEMENT IS FORMALLY
C             CONSIDERED AS TWO SIMPLE ELEMENTS.
C
C             OFTEN IT IS CONVENIENT TO WORK WITH WAVES, THE RAYS OF
C             WHICH HAVE EITHER A COMPOUND ELEMENT IN A COMPLEX BLOCK OR
C             TWO SIMPLE ELEMENTS OF THE SAME WAVE TYPE (UNCONVERTED
C             REFLECTION) IN THE SAME BLOCK, AS WITH ONE ELEMENTARY
C             WAVE.  THE COMPOUND-ELEMENT BLOCK CODE MAKES A DIVISION OF
C             THE WAVE FIELD INTO SUCH ELEMENTARY WAVES POSSIBLE.
C
C             WE INTRODUCE THE COMPOUND-ELEMENT BLOCK CODE AS FOLLOWS.
C             FOR ANY SIMPLE ELEMENT OF A RAY ONE NUMBER IS USED.  THE
C             ABSOLUTE VALUE OF THIS NUMBER SPECIFIES THE INDEX OF THE
C             COMPLEX BLOCK, IN WHICH THE SIMPLE ELEMENT IS SITUATED.
C             THE SIGN OF THIS NUMBER SPECIFIES THE TYPE OF THE WAVE
C             (PLUS SIGN...  P WAVE, MINUS SIGN...  S WAVE).  FOR ANY
C             COMPOUND ELEMENT OF A RAY, TWO IDENTICAL NUMBERS ARE USED.
C             A COMPOUND ELEMENT AND A DOUBLET OF SIMPLE ELEMENTS,
C             DESCRIBED BY THE SAME NUMBERS, ARE NOT DISTINGUISHED.  THE
C             CODE OF A RAY IS A CHAIN OF THESE NUMBERS DESCRIBING THE
C             RAY FROM ITS INITIAL POINT TO ITS ENDPOINT.  THE USE OF
C             THE COMPOUND-ELEMENT BLOCK CODE LEADS TO MORE GENERAL
C             ELEMENTARY WAVES.
C
C     (4) GENERALIZED LAYER CODE
C             ANOTHER CODE MAY BE OBTAINED FROM THE COMPOUND-ELEMENT
C             BLOCK CODE IF THE MODIFIED INTERPRETATION OF THE MODEL
C             STRUCTURE IS USED.  THE MODIFIED INTERPRETATION CONSISTS
C             IN ASSUMING AN EXISTENCE OF FICTITIOUS PARTS OF BLOCKS OF
C             A ZERO THICKNESS SITUATED BETWEEN EVERY NEIGHBOURING
C             COMPLEX BLOCKS, THE INDICES OF WHICH ARE NOT SEQUENTIALLY
C             ORDERED.  THE NUMBER OF FICTITIOUS BLOCKS IS CONSIDERED TO
C             BE JUST THE NUMBER, WHICH IS NECESSARY TO FILL THE GAP
C             BETWEEN INDICES OF THE TWO NEIGHBOURING COMPLEX BLOCKS.
C             THE RAY CAN ONLY PASS THROUGH THE FICTITIOUS BLOCKS, NO
C             REFLECTION IN FICTITIOUS BLOCKS BEING ALLOWED.  SIMILARLY,
C             NO CONVERSION IS ALLOWED NEITHER AT INTERFACES BETWEEN
C             FICTITIOUS BLOCKS NOR AT THE INTERFACE AT WHICH THE RAY
C             LEAVES A FICTITIOUS BLOCK (CONVERSION MAY TAKE PLACE ONLY
C             AT THE INTERFACE AT WHICH THE RAY ENTERS THE FICTITIOUS
C             BLOCK).  IF ANY OF THESE PROHIBITED SITUATIONS IS
C             SPECIFIED BY THE CODE, IT LEADS TO THE TERMINATION OF
C             COMPUTATIONS OF THE CORRESPONDING RAY.
C
C             THIS CODE IS A GENERALIZATION OF THE CODE USED IN THE 2-D
C             PROGRAM PACKAGE 'SEIS83' AND IS DESCRIBED IN CERVENY,
C             MOLOTKOV AND PSENCIK: RAY METHOD IN SEISMOLOGY, CHARLES
C             UNIVERSITY, PRAGUE 1977.  IT ENABLES ANY BLOCK STRUCTURE
C             TO BE INTERPRETED AS LOCALLY LAYERED STRUCTURE WITH
C             FICTITIOUS LAYERS.  THUS, THE ROUTINES FOR AUTOMATIC OR
C             SEMIAUTOMATIC GENERATION OF RAY CODES FOR LAYERED
C             STRUCTURES MAY BE USED FOR THIS CODE EVEN IN GENERAL BLOCK
C             STRUCTURES.  ITS EFFECTIVE USE IS CONDITIONED BY PROPER
C             INDEXING OF COMPLEX BLOCKS, WHICH MINIMIZES THE NUMBER OF
C             FICTITIOUS LAYERS.  FOR EXAMPLE, IN A LAYERED STRUCTURE,
C             THE BLOCKS (LAYERS) SHOULD BE INDEXED SEQUENTIALLY FROM
C             THE TOP TO THE BOTTOM.
C
C     (5) REFLECTION-TRANSMISSION CODE
C             FOR ANY ELEMENT OF A RAY ONE INTEGER IS USED.  THE FIRST
C             ELEMENT OF THE RAY, I.E.  THE ELEMENT CONTAINING THE
C             INITIAL POINT OF THE RAY, IS DENOTED BY 1 FOR P AND BY -1
C             FOR S WAVE.  ANY OTHER ELEMENT IS DENOTED BY 1 OR -1 (P OR
C             S WAVE) IF THE RAY IS TRANSMITTED AT THE INITIAL POINT OF
C             THE ELEMENT, AND BY 2 OR -2 (P OR S WAVE) IF THE RAY IS
C             REFLECTED AT THE INITIAL POINT OF THE ELEMENT.  THE CODE
C             OF A RAY IS A CHAIN OF THE ABOVE NUMBERS CORRESPONDING TO
C             INDIVIDUAL ELEMENTS OF THE RAY FROM ITS INITIAL POINT TO
C             ITS ENDPOINT.  IN A LAYERED MODEL, THIS CODE IS VERY
C             CONVENIENT FOR THE DESCRIPTION OF REFRACTED AND PRIMARY
C             REFLECTED WAVES.  THE APPLICATION OF THIS CODE IS
C             STRAIGHTFORWARD.
C
C SPECIFICATION OF AN ELEMENTARY WAVE
C     AN ELEMENTARY WAVE IS SPECIFIED BY THE FOLLOWING DATA:
C     (1) INTEGER KODTYP, WHICH DETERMINES THE TYPE OF THE CODE
C             1... BLOCK-SURFACE CODE,
C             2... BASIC BLOCK CODE,
C             3... COMPOUND-ELEMENT BLOCK CODE,
C             4... GENERALIZED LAYER CODE,
C             5... TRANSMISSION-REFLECTION CODE.
C     (2) INTEGER ARRAY KODE CONTAINING THE CODE OF THE ELEMENTARY
C             WAVE.  THE CODE IS A FINITE SEQUENCE OF NONZERO INTEGERS.
C             A ZERO INDICATES THE END OF A CODE IN THE INPUT DATA.
C     THE DATA (1) AND (2) SHOULD BE STORED IN COMMON BLOCK
C     --------------------------------
C     COMMON /COD/ KODTYP,KODE(MKODE).
C     --------------------------------
C     THE DIMENSION MKODE OF THE ARRAY KODE MAY BE ADJUSTED BY THE USER.
C
C DATE: 1988, JUNE 3
C WRITTEN BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
C
C=======================================================================
C
      SUBROUTINE CODE1(LUN,IWAVE,IWAVE0,IKODE)
      INTEGER LUN,IWAVE,IWAVE0,IKODE
C
C THIS SUBROUTINE IS CALLED WHEN STARTING THE COMPUTATION OF A NEW
C ELEMENTARY WAVE.  IT STORES THE CODE OF THE ELEMENTARY WAVE INTO THE
C COMMON BLOCK /COD/.
C
C INPUT:
C     LUN...  LOGICAL UNIT NUMBER OF THE EXTERNAL INPUT DEVICE
C             CONTAINING THE INPUT DATA.
C     IWAVE...ZERO WHEN STARTING THE COMPLETE RAY TRACING PROGRAM,
C             OTHERWISE THE INDEX OF THE LAST ELEMENTARY WAVE WHICH HAS
C             BEEN COMPUTED (I.E. THE OUTPUT FROM THE PREVIOUS
C             INVOCATION OF THE SUBROUTINE CODE1).
C
C OUTPUT:
C     IWAVE...ZERO IF ALL REQUIRED ELEMENTARY WAVES ARE COMPUTED AND THE
C             COMPLETE RAY TRACING PROGRAM WILL BE TERMINATED.
C             OTHERWISE, THE INDEX OF THE ELEMENTARY WAVE WHICH WILL BE
C             COMPUTED (I.E. NSKIP+1 FOR THE FIRST COMPUTED ELEMENTARY
C             WAVE, WHERE NSKIP IS THE NUMBER OF ELEMENTARY WAVES HAVING
C             BEEN SKIPPED OVER, OTHERWISE THE INPUT VALUE INCREASED BY
C             ONE).
C     IWAVE0..INDEX OF THE ALREADY COMPUTED ELEMENTARY WAVE HAVING THE
C             MOST NUMEROUS COMMON ELEMENTS WITH THE CURRENT ELEMENTARY
C             WAVE.  IN THE CASE OF SEVERAL POSSIBILITIES, THE FIRST
C             COMPUTED WAVE OF THEM IS TAKEN.
C     IKODE...THE LENGTH OF THE COMMON PART OF THE CODES OF THE IWAVE-TH
C             AND IWAVE0-TH ELEMENTARY WAVES.
C
C COMMON BLOCK /COD/:
      INTEGER MKODE
      PARAMETER (MKODE=128)
      INTEGER KODTYP,KODE(MKODE)
      COMMON /COD/ KODTYP,KODE
C ALL THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE DEFINED IN THIS
C SUBROUTINE.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL CODEB,LCODE
C     CODEB.. BLOCK DATA SUBROUTINE OF THIS FILE.
C     LCODE.. THIS FILE.
C
C ERROR MESSAGES:
C     401...  INSUFFICIENT MEMORY FOR CODES:
C             THE DIMENSION MKODES OF THE ARRAY KODES(0:MKODES) IN
C             THIS SUBROUTINE SHOULD BE INCREASED.  THE DIMENSION MKODES
C             SHOULD AT LEAST EQUAL TO THE NUMBER OF ELEMENTARY WAVES
C             PLUS THE TOTAL NUMBER OF ALL INDECES FORMING THE CODES OF
C             THE WAVES.
C             NOTE:
C             IF THE PARAMETER MODCRT OF THE DATA SET 'DCRT' (SEE SOURCE
C             CODE FILE 'RAY.FOR') IS ZERO, THE OUTPUT PARAMETERS IWAVE0
C             AND IKODE OF THIS SUBROUTINE ARE LIKELY UNUSED.  THEN,
C             INSTEAD OF INCREASING MKODES, THE PAUSE AND STOP
C             STATEMENTS GENERATING THIS ERROR MESSAGE MAY BE DISABLED
C             TO SAVE THE COMPUTER MEMORY.  THEN THIS SUBROUTINE WILL
C             WORK, BUT CONSIDERING ONLY THE FIRST SEVERAL ELEMENTARY
C             WAVES WHEN CALCULATING OUTPUT VALUES OF IWAVE0 AND IKODE.
C
C DATE: 1993, DECEMBER 18
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     STORAGE FOR PREVIOUS CODES:
      INTEGER MKODES
      PARAMETER (MKODES=1024)
      INTEGER KODES(0:MKODES)
      SAVE KODES
C
C     AUXILIARY STORAGE LOCATIONS:
      CHARACTER*80 TEXTC
      INTEGER NSKIP,I,J,K,L
      SAVE NSKIP
C
C     TEXTC...THE NAME OF THE DATA. STRING OF 80 CHARACTERS.
C     NSKIP...THE NUMBER OF ELEMENTARY WAVES TO BE SKIPPED OVER.
C     I,J...  AUXILIARY LOOP VARIABLES.
C     K...    AUXILIARY VARIABLE - INDEX IN ARRAY KODES.
C     L...    AUXILIARY VARIABLE - LENGTH OF THE CURRENT CODE.
C
C.......................................................................
C
      IF(IWAVE.EQ.0) THEN
C       READING THE NAME OF THE INPUT DATA
        READ(LUN,*) TEXTC
C
C       READING THE NUMBER OF ELEMENTARY WAVES TO BE SKIPPED OVER
        NSKIP=0
        READ(LUN,*) NSKIP
C       NO CODES OF ELEMENTARY WAVES ARE READ AND STORED NOW
        KODES(0)=0
      END IF
C
C     READING THE CODE OF THE ELEMENTARY WAVE
   10 CONTINUE
        KODTYP=0
        DO 11 I=1,MKODE
          KODE(I)=0
   11   CONTINUE
        READ(LUN,*) KODTYP,KODE
        IF(KODTYP.EQ.0) THEN
          IWAVE=0
          RETURN
        ELSE
          K=KODES(IWAVE)+1
          L=LCODE()
          IWAVE=IWAVE+1
          IF(K+L.GT.MKODES) THEN
            PAUSE 'ERROR 401 IN CODE1: INSUFFICIENT MEMORY FOR CODES'
            STOP
          ELSE
C           STORING THE CODES
            DO 12 I=L,1,-1
              KODES(K+I)=KODE(I)
   12       CONTINUE
            DO 13 I=K-1,IWAVE,-1
              KODES(I+1)=KODES(I)
   13       CONTINUE
            KODES(IWAVE)=K+L
            DO 14 I=IWAVE-1,0,-1
              KODES(I)=KODES(I)+1
   14       CONTINUE
          END IF
        END IF
      IF(IWAVE.LE.NSKIP) GO TO 10
C
C     DETERMINING ELEMENTARY WAVE HAVING THE MOST NUMEROUS COMMON
C     ELEMENTS WITH THE CURRENT ELEMENTARY WAVE
      IKODE=0
      IWAVE0=0
      L=LCODE()
      DO 39 J=1,MIN0(KODES(0),IWAVE-1)
        K=KODES(J-1)
        DO 31 I=1,MIN0(KODES(J)-K,L)
          IF(KODE(I).NE.KODES(K+I)) THEN
            GO TO 32
          END IF
   31   CONTINUE
   32   CONTINUE
        I=I-1
        IF(I.GT.IKODE) THEN
          IKODE=I
          IWAVE0=J
        END IF
   39 CONTINUE
      RETURN
      END
C
C=======================================================================
C
      INTEGER FUNCTION LCODE()
C
C INTEGER FUNCTION KOOR IS DESIGNED TO RETURN THE LENGTH OF THE CODE OF
C THE CURRENT ELEMENTARY WAVE.
C
C NO INPUT.
C
C OUTPUT:
C     LCODE...LENGTH OF THE CODE OF THE CURRENT ELEMENTARY WAVE, I.E.
C             THE COUNT OF THE NUMERIC ITEMS IN THE ARRAY KODE WHICH
C             DESCRIBES THE BEHAVIOUR OF A RAY AT INTERFACES.  THE ARRAY
C             IS CALLED HERE THE CODE OF ELEMENTARY WAVE.
C
C COMMON BLOCK /COD/:
      INTEGER MKODE
      PARAMETER (MKODE=128)
      INTEGER KODTYP,KODE(MKODE)
      COMMON /COD/ KODTYP,KODE
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 I
C
      DO 1 I=1,MKODE
        IF(KODE(I).EQ.0) THEN
          LCODE=I-1
          GO TO 2
        END IF
    1 CONTINUE
      LCODE=MKODE
    2 CONTINUE
      RETURN
      END
C
C=======================================================================
C
      INTEGER FUNCTION NELEM(KODIND)
      INTEGER KODIND
C
C THIS FUNCTION IS DESIGNED TO RETURN THE MAXIMUM POSSIBLE NUMBER OF RAY
C ELEMENTS CORRESPONDING TO THE GIVEN POSITION IN THE CODE OF THE
C ELEMENTARY WAVE.
C
C INPUT:
C     KODIND..POSITION IN THE CODE (INDEX IN ARRAY KODE)
C NONE OF THE INPUT PARAMETERS ARE ALTERED.
C
C OUTPUT:
C     NELEM...NUMBER OF RAY ELEMENTS BETWEEN THE INITIAL POINT OF THE
C             RAY AND THE END OF THE ELEMENT CORRESPONDING TO THE GIVEN
C             POSITION IN THE CODE.  NOTE THAT, FOR A PARTICULAR RAY,
C             SOME OF THESE POSSIBLE RAY ELEMENTS MAY BE LEFT OUT OR
C             COUPLED INTO ONE VIRTUAL ELEMENT.
C
C COMMON BLOCK /COD/:
      INTEGER MKODE
      PARAMETER (MKODE=128)
      INTEGER KODTYP,KODE(MKODE)
      COMMON /COD/ KODTYP,KODE
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
      IF(KODTYP.EQ.1) THEN
        NELEM=(KODIND+1)/2
      ELSE
        NELEM=KODIND
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE CODE(IY,KODIND,ICBNEW,IEND)
      INTEGER IY(12),KODIND,ICBNEW,IEND
C
C THIS SUBROUTINE TRANSFORMS THE USED NUMERICAL CODE OF THE ELEMENTARY
C WAVE UNDER CONSIDERATION INTO INSTRUCTIONS SPECIFYING THE BEHAVIOUR
C OF THE WAVE AT THE INITIAL POINTS OF RAYS AND AT ALL POINTS OF
C INCIDENCE OF THE RAYS AT INTERFACES (BOUNDARIES OF COMPLEX BLOCKS).
C THUS, IT IS FIRST CALLED AT THE INITIAL POINT OF A RAY AND THEN
C SUCCESSIVELY AT ALL POINTS OF INCIDENCE.
C
C INPUT:
C     IY...   INTEGER ARRAY. ITS NUMERICAL STORAGE UNITS 2, 3, 5, 6, 8
C             MUST BE DEFINED AS FOLLOWS:
C             IY(2)... POSITION IN THE CODE (INDEX IN ARRAY KODE)
C               CORRESPONDING TO THE LAST COMPUTED ELEMENT OF A RAY.
C               IY(2)=0 AT THE INITIAL POINT OF THE RAY,
C               IY(2)=KODIND FROM THE LAST INVOCATION OF SUBROUTINE CODE
C                 AT OTHER POINTS OF THE RAY.
C             IY(3)=ICB0... INDEX OF THE NEIGHBOURING COMPLEX BLOCK,
C               FROM WHICH THE RAY ENTERED THE COMPLEX BLOCK IN WHICH
C               THE LAST COMPUTED ELEMENT OF THE RAY IS SITUATED.
C               IY(3)=0 BEFORE LEAVING THE COMPLEX BLOCK, IN WHICH THE
C               INITIAL POINT OF THE RAY IS SITUATED.
C               AT THE INITIAL POINT OF THE RAY (IY(2)=0), IY(3) IS
C               IGNORED.
C             IY(5)=ICB1... INDEX OF THE COMPLEX BLOCK CONTAINING THE
C               COMPUTED ELEMENT OF THE RAY, SUPPLEMENTED BY THE SIGN
C               '+' FOR P WAVE AND SIGN '-' FOR S WAVE.
C               ICB1 IS IGNORED AT THE INITIAL POINT OF THE RAY
C               (IY(2)=0).
C             IY(6)=ISRF... INDEX OF THE SURFACE AT WHICH THE ENDPOINT
C               OF THE LAST COMPUTED ELEMENT OF THE RAY IS SITUATED.
C               THE SIGN OF IY(6) IS IGNORED.
C               IY(6)=0 AT THE INITIAL POINT OF THE RAY.
C             IY(8)=ICB2... INDEX OF THE COMPLEX BLOCK TOUCHING THE
C               COMPLEX BLOCK ICB1 FROM THE OTHER SIDE OF THE SURFACE
C               ISRF AT THE ENDPOINT OF THE LAST COMPUTED ELEMENT OF THE
C               RAY. IY(8)=0 FOR A FREE SPACE.
C               AT THE INITIAL POINT OF THE RAY, ICB2 IS THE INDEX OF
C               THE COMPLEX BLOCK CONTAINING THE INITIAL POINT.
C THE INPUT PARAMETER IS NOT ALTERED.
C
C OUTPUT:
C     KODIND..POSITION IN THE CODE (INDEX IN THE ARRAY KODE) CORRES-
C             PONDING TO THE NEXT ELEMENT OF THE RAY.
C     ICBNEW..THE INDEX OF THE COMPLEX BLOCK IN WHICH THE NEXT
C             ELEMENT OF THE RAY IS TO BE SITUATED, SUPPLEMENTED
C             BY THE SIGN "+" FOR P WAVE OR "-" FOR S WAVE.
C     IEND... INFORMATION ON THE PROCESS OF THE INTERPRETATION
C             OF THE CODE:
C             IEND.EQ.0... COMPUTATION OF THE RAY CONTINUES.
C             IEND.NE.0... THE COMPUTATION OF THE RAY TERMINATES.
C               DIFFERENT VALUES OF IEND SPECIFY THE REASON FOR THE
C               TERMINATION:
C               10... RAY SATISFIES THE WHOLE CODE (REGULAR TERMINATION
C                 OF THE RAY).
C               21... THE POINT OF INCIDENCE IS SITUATED AT A DIFFERENT
C                 SURFACE THAN THAT REQUIRED BY THE CODE.
C               22... THE NEXT ELEMENT OF THE RAY IS REQUIRED BY THE
C                 CODE TO BE SITUATED IN A COMPLEX BLOCK THAT DOES NOT
C                 TOUCH THE POINT OF INCIDENCE.
C               23... TRANSMISSION IS REQUIRED BY THE CODE AT A FREE
C                 SURFACE.
C               30... REFLECTION OR WAVE CONVERSION AT THE FICTITIOUS
C                 PART OF THE INTERFACE.
C
C COMMON BLOCK /COD/:
C     ------------------------------------------------------------------
      INTEGER MKODE
      PARAMETER (MKODE=128)
      INTEGER KODTYP,KODE(MKODE)
      COMMON /COD/ KODTYP,KODE
C     ------------------------------------------------------------------
C     KODTYP..DETERMINES THE TYPE OF THE CODE:
C             1... BLOCK-SURFACE CODE,
C             2... BASIC BLOCK CODE,
C             3... COMPOUND-ELEMENT BLOCK CODE,
C             4... GENERALIZED LAYER CODE,
C             5... TRANSMISSION-REFLECTION CODE.
C     KODE... ARRAY CONTAINING THE CODE OF THE ELEMENTARY WAVE.
C             THE CODE IS A FINITE SEQUENCE OF NONZERO INTEGERS.
C             ZERO INDICATES THE END OF A CODE.
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED.
C
C DATE: 1990, NOVEMBER 18
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER I,J
C
      KODIND=IY(2)
      GO TO (10,20,30,30,50),KODTYP
C
C     BLOCK-SURFACE CODE:
   10 CONTINUE
      IF(KODIND.GT.0) THEN
C       CHECK ON THE SURFACE:
        KODIND=KODIND+1
        IF(KODE(KODIND).EQ.0) THEN
C         END OF CODE
          IEND=10
          RETURN
        ELSE IF(KODE(KODIND).NE.IABS(IY(6))) THEN
C         WRONG SURFACE
          IEND=21
          RETURN
        END IF
      END IF
C     THE REST OF INTERPRETATION IS THE SAME AS IN THE BASIC BLOCK CODE
C
C     BASIC BLOCK CODE:
   20 CONTINUE
      KODIND=KODIND+1
      ICBNEW=KODE(KODIND)
      IF(ICBNEW.EQ.0) THEN
C       END OF CODE
        IEND=10
      ELSE IF(IABS(ICBNEW).EQ.IY(8)) THEN
C       TRANSMISSION
        IEND=0
      ELSE IF(KODIND.GT.1.AND.IABS(ICBNEW).EQ.IABS(IY(5))) THEN
C       REFLECTION
        IEND=0
      ELSE
C       REQUIRED COMPLEX BLOCK IS NOT ATTAINABLE
        IEND=22
      END IF
      RETURN
C
C     COMPOUND-ELEMENT BLOCK CODE:
   30 CONTINUE
      IF(KODIND.GT.0) THEN
        I=IABS(IY(5))
        J=(I-IABS(IY(3)))*(I-IY(8))
        DO 31 I=KODIND-1,1,-1
          IF(KODE(I).NE.KODE(KODIND)) GO TO 32
          J=-J
   31   CONTINUE
   32   CONTINUE
        IF(J.GT.0) THEN
C         COMPOUND ELEMENT:
          KODIND=KODIND+1
          IF(KODE(KODIND).NE.IY(5)) THEN
C           WRONG SURFACE
            IEND=21
            RETURN
          END IF
        END IF
      END IF
C     REST OF INTERPRETATION IS THE SAME AS IN BASIC BLOCK CODE
      IF(KODTYP.EQ.3) GO TO 20
C
C     GENERALIZED LAYER CODE:
C     INTERPRETATION OF COMPOUND ELEMENTS IS THE SAME AS IN
C     COMPOUND-ELEMENT BLOCK CODE, THEN:
      KODIND=KODIND+1
      ICBNEW=KODE(KODIND)
      IF(ICBNEW.EQ.0) THEN
C       END OF CODE
        IEND=10
      ELSE IF(KODIND.EQ.1) THEN
C       INITIAL POINT OF A RAY:
        IF(IABS(ICBNEW).EQ.IY(8)) THEN
C         INITIAL POINT OF A RAY IS SITUATED IN THE REQUIRED C.B.
          IEND=0
        ELSE
C         INITIAL POINT OF A RAY IS NOT SITUATED IN THE REQUIRED C.B.
          IEND=22
        END IF
      ELSE IF(IABS(ICBNEW).EQ.IABS(IY(5))) THEN
C       REFLECTION
        IEND=0
      ELSE
C       POSSIBLE TRANSMISSION:
        J=ISIGN(1,IY(8)-IABS(IY(5)))
        IF(IABS(ICBNEW).EQ.IABS(IY(5))+J) THEN
C         LOOP FOR FICTITIOUS PARTS OF BLOCKS:
          DO 41 I=IABS(IY(5))+J+J,IY(8),J
            KODIND=KODIND+1
            ICBNEW=KODE(KODIND)
            IF(ICBNEW.EQ.0) THEN
C             END OF CODE IN THE FICTITIOUS PART OF BLOCK
              IEND=10
              RETURN
            ELSE IF(ICBNEW.EQ.ISIGN(I,KODE(KODIND-1))) THEN
C             TRANSMISSION FROM THE FICTITIOUS PART OF BLOCK
              IEND=0
            ELSE
C             TERMINATION OF THE RAY COMPUTATION IN THE FICTITIOUS PART
C             OF THE BLOCK
              IEND=30
              RETURN
            END IF
   41     CONTINUE
        ELSE
C         REQUIRED COMPLEX BLOCK IS NOT ATTAINABLE
          IEND=22
        END IF
      END IF
      RETURN
C
C     TRAMSMISSION-REFLECTION CODE:
   50 CONTINUE
      KODIND=KODIND+1
      I=KODE(KODIND)
      IF(IABS(I).EQ.1) THEN
C       TRANSMISSION:
        IF(IY(8).GT.0) THEN
C         TRANSMISSION INTO MATERIAL BLOCK
          ICBNEW=ISIGN(IY(8),I)
          IEND=0
        ELSE
C         TRANSMISSION INTO FREE SPACE
          IEND=23
        END IF
      ELSE IF(IABS(I).EQ.2.AND.KODIND.GT.1) THEN
C       REFLECTION
        ICBNEW=ISIGN(IABS(IY(5)),I)
        IEND=0
      ELSE
C       END OF CODE
        IEND=10
      END IF
      RETURN
      END
C
C=======================================================================
C
