C PROGRAM 'GRID' TO GENERATE THE VELOCITIES IN A RECTANGULAR GRID
C REQUIRED FOR THE FULL WAVE FINITE DIFFERENCES, THE SHORTEST PATH
C CALCULATION OF SEISMIC RAYS, EIKONAL EQUATION 'FINITE DIFFERENCES',
C OR RASTER IMAGING OF THE MODEL.  ALSO AN OBLIQUE VERTICAL 2-D SECTION
C ACROSS THE 3-D MODEL MAY BE GRIDDED.
C
C THE RECTANGULAR GRID IS SPECIFIED IN CARTESIAN COORDINATES, AND THEN
C TRANSFORMED TO THE MODEL COORDINATES.  FOR THE CARTESIAN COORDINATES
C CONNECTED WITH A PARTICULAR KIND OF CURVILINEAR MODEL COORDINATES SEE
C SUBROUTINE CARTES OF THE FILE 'METRIC.FOR'.
C THE RECTANGULAR GRID COULD ALSO BE SPECIFIED IN RESPECT TO THE MODEL
C COORDINATES, AND LIMITED BY COORDINATE PLANES SPECIFIED IN THE MODEL
C COORDINATES.  THIS OPTION MAY BE ENABLED BY CHANGING THE FIRST
C EXECUTIVE STATEMENT '      LCART=.TRUE.' TO '      LCART=.FALSE.'.
C
C NAME OF THE MAIN INPUT DATA FILE READ FROM THE * EXTERNAL UNIT:
C (1) 'GRID'
C     'GRID'..NAME OF THE MAIN INPUT DATA FILE DESCRIBED BELOW.
C     DEFAULT:  'GRID'='GRID.DAT'.
C
C MAIN INPUT DATA FILE 'GRID':
C (1) 'MODEL','NET','IND','VEL','ICB','VELDER',/
C     'MODEL'... INPUT DATA FILE DESCRIBING THE MODEL, IT IS DESCRIBED
C             IN THE SUBROUTINE FILE 'MODEL.FOR'.
C     'NET'...NAME OF THE INPUT DATA FILE, DESCRIBED BELOW, DESCRIBING
C             THE DIMENSIONS AND DENSITY OF A RECTANGULAR GRID OF
C             POINTS.
C     'IND'...IF NOT BLANK, THE NAME OF THE INDEX FILE.
C             THIS OPTION ENABLES TO SPECIFY OTHER THAN RECTANGULAR
C             REGION COVERED BY A RECTANGULAR GRID:
C             THE RECTANGULAR VOLUME BOUNDED BY COORDINATE LIMITS
C             X1MIN,X1MAX, X2MIN,X2MAX, AND X3MIN,X3MAX IS DIVIDED INTO
C             N1*N2*N3 BIG BRICKS.  EACH ELEMENT (INDEX) OF THE INDEX
C             FILE CORRESPONDS TO ONE BIG BRICK.  IF IT EQUALS ZERO,
C             THE BIG BRICK DOES NOT BELONG TO THE REGION IN WHICH THE
C             VELOCITY IS DISCRETIZED.
C     'VEL'...NAME OF THE OUTPUT FORMATTED FILE,
C             CONTAINING THE VELOCITIES AT THE GIVEN GRIDPOINTS.
C             VELOCITY=0 IS INSERTED IN A FREE SPACE.
C             IF BLANK, THE FILE IS NOT CREATED.
C     'ICB'...NAME OF THE OUTPUT FORMATTED FILE,
C             CONTAINING THE INDICES OF COMPLEX GEOLOGICAL BLOCKS AT THE
C             GIVEN GRIDPOINTS.
C             IF BLANK, THE FILE IS NOT CREATED.
C     'VELDER'... NAME OF THE OUTPUT FORMATTED FILE,
C             CONTAINING VELOCITIES WITH THEIR FIRST AND SECOND
C             DERIVATIVES AT THE GIVEN GRIDPOINTS.
C             IF BLANK, THE FILE IS NOT CREATED.
C     DEFAULT:  'MODEL'='MODEL.DAT', 'NET'='NET.DAT', 'IND'=' ',
C             'VEL'='VEL.OUT', 'ICB'='ICB.OUT', 'VELDER'=' '.
C (2) (IPS(I),I=1,NPS),/
C     IPS...  LIST OF INDICES OF COMPLEX BLOCKS TERMINATED BY A SLASH.
C             IF A COMPLEX BLOCK IS NOT LISTED, P-WAVE VELOCITY IS
C             EVALUATED AT GRID POINTS.
C             IPS(I).LT.0:  S-WAVE VELOCITY IS EVALUATED AT GRID POINTS.
C             IPS(I).GT.0:  FREE SPACE (VELOCITY=0) IS ASSUMED AT GRID
C               POINTS.  THIS OPTION IS LIKELY BE USED TO AVOID
C               REFRACTION IN DEEP LAYERS WHEN COMPUTING REFLECTED RAYS
C               BY MEANS OF THE SHORTEST PATH NETWORK ALGORITHM.
C             DEFAULT: P-WAVE VELOCITY IN ALL MATERIAL COMPLEX BLOCKS.
C
C INPUT FILE 'MODEL' IS DESCRIBED WITHIN THE FORTRAN77 SOURCE CODE FILE
C     'MODEL.FOR'.
C
C INPUT FILE 'NET':
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, THE INPUT PARAMETER IS OF THE TYPE REAL.
C (1) ONE OF THE FOLLOWING POSSIBILITIES (1.1), (1.2), (1.3), OR (1.4):
C (1.1) N1,N2,N3,/
C             A BASIC OPTION:
C     N1,N2,N3... NUMBERS OF GRIDPOINTS OF THE RECTANGULAR GRID IN THE
C             DIRECTIONS OF AXES X1,X2,X3.
C             THE RECTANGULAR VOLUME BOUNDED BY COORDINATE LIMITS X1MIN,X1MAX,
C             X2MIN,X2MAX, AND X3MIN,X3MAX IS DIVIDED INTO N1*N2*N3 BIG BRICKS.
C             IF THE NUMBERS L1,L2,L3 ARE SPECIFIED IN ADDITION TO N1,N2,N3,
C                     EACH BIG BRICK IS SUBDIVIDED INTO L1*L2*L3 SMALL BRICKS.
C             IF THE NUMBERS L1,L2,L3 ARE NOT SPECIFIED, EACH BIG BRICK CONTAINS
C                     JUST ONE SMALL BRICK, AS LARGE AS BIG ONE.
C             GRIDPOINTS ARE DEFINED AS CENTRES OF THE SMALL BRICKS.
C     DEFAULT: N3=1.
C (1.2) N1,N2,N3,L1,L2,L3,/
C             A SPECIAL OPTION:
C             THE RECTANGULAR VOLUME BOUNDED BY COORDINATE LIMITS
C             X1MIN,X1MAX, X2MIN,X2MAX, AND X3MIN,X3MAX IS DIVIDED INTO
C             N1*N2*N3 BIG BRICKS.
C             THEN EACH BIG BRICK IS DIVIDED INTO L1*L2*L3 SMALL BRICKS.
C             THE OUTPUT VELOCITIES ARE COMPUTED IN THE CENTRES OF SMALL
C             BRICKS.
C             OUTER LOOP IS OVER BIG BRICKS, THE DISCRETE VELOCITIES
C             WITHIN EACH BIG BRICK BEING OUTPUT CONSECUTIVELY.
C     N1,N2,N3... NUMBERS OF BIG BRICKS IN THE X1,X2,X3-AXIS DIRECTIONS.
C             MUST BE POSITIVE.
C             N3 NEED NOT BE SPECIFIED FOR A 2-D MODEL.
C     L1,L2,L3... NUMBERS OF SMALL BRICKS IN ONE BIG BRICK IN THE
C             X1,X2,X3-AXIS DIRECTIONS.  IF SPECIFIED, MUST BE POSITIVE.
C     DEFAULT: N3=1, L1=1, L2=1, L3=1.
C (1.3) 0,N2,N3,/
C             A SPECIAL OPTION:  2-D OBLIQUE VERTICAL SECTION IN 3-D:
C             THE RECTANGULAR VERTICAL SECTION BOUNDED BY THE VERTICAL
C             LINES (X1,X2)=(X1MIN,X2MIN) AND (X1,X2)=(X1MAX,X2MAX),
C             FROM X3=X3MIN TO X3=X3MAX IS DIVIDED INTO N2*N3 CELLS.
C     N2,N3.. NUMBERS OF CELLS IN THE HORIZONTAL AND VERTICAL
C             DIRECTIONS.  MUST BE POSITIVE.
C     DEFAULT: N3=1.
C (1.4) 0,N2,N3,L1,L2,L3,/
C             A SPECIAL OPTION:  2-D OBLIQUE VERTICAL SECTION IN 3-D:
C             THE RECTANGULAR VERTICAL SECTION BOUNDED BY THE VERTICAL
C             LINES (X1,X2)=(X1MIN,X2MIN) AND (X1,X2)=(X1MAX,X2MAX),
C             FROM X3=X3MIN TO X3=X3MAX IS DIVIDED INTO N2*N3 BIG CELLS,
C             EACH BIG CELL BEING DIVIDED INTO L2*L3 SMALL CELLS.
C     N2,N3.. NUMBERS OF BIG CELLS IN THE HORIZONTAL AND VERTICAL
C             DIRECTIONS.  MUST BE POSITIVE.
C     L1...   MUST EQUAL 0 OR 1.  L1=0 IS UNDERSTOOD AS L1=1.
C     L2,L3.. NUMBERS OF SMALL CELLS IN ONE BIG CELL IN THE HORIZONTAL
C             AND VERTICAL DIRECTIONS.  IF SPECIFIED, MUST BE POSITIVE.
C     DEFAULT: N3=1, L1=1, L2=1, L3=1.
C (2) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,/
C     BOUNDARIES OF THE VOLUME COVERED BY THE RECTANGULAR HOMOGENEOUS
C     GRID.  THE X1 COORDINATES OF THE GRIDPOINTS ARE
C             X1MIN+0.5*DX1, X1MIN+1.5*DX1, ... ,X1MAX-0.5*DX1  ,
C     WHERE
C             DX1=(X1MAX-X1MIN)/N1,
C     OR
C             DX1=(X1MAX-X1MIN)/N1*L1,
C     IF L1 IS SPECIFIED, SIMILARLY FOR X2 AND X3.
C     DEFAULT: X1MIN=0, X1MAX=0, X2MIN=0, X2MAX=0, X3MIN=0, X3MAX=0.
C
C INPUT FILE 'IND':
C             THIS OPTION ENABLES TO SPECIFY OTHER THAN RECTANGULAR
C             REGION COVERED BY A RECTANGULAR GRID:
C             THE RECTANGULAR VOLUME BOUNDED BY COORDINATE LIMITS
C             X1MIN,X1MAX, X2MIN,X2MAX, AND X3MIN,X3MAX IS DIVIDED INTO
C             N1*N2*N3 BIG BRICKS.  EACH ELEMENT (INDEX) OF THE INDEX
C             FILE CORRESPONDS TO ONE BIG BRICK.  IF IT EQUALS ZERO,
C             THE BIG BRICK DOES NOT BELONG TO THE REGION IN WHICH THE
C             VELOCITY IS DISCRETIZED.
C     ATTENTION: THE NONZERO INDICES MUST BE FORMED BY THE SEQUENCE
C             1,2,3,... OF POSITIVE INTEGERS.
C (1) (IND(I),I=1,N1*N2*N3)
C     IND(I)..ZERO IF THE I-TH BIG BRICK DOES NOT BELONG TO THE REGION
C               IN WHICH THE VELOCITY IS DISCRETIZED.
C             OTHERWISE THE INDEX THE BIG BRICK.  THE GRIDPOINTS WITHIN
C               THE BIG BRICK ARE INDEXED BY INTEGERS FROM
C               L1*L2*L3*(IND(I)-1)+1 TO L1*L2*L3*IND(I).
C     DEFAULT: IND(I)=I.
C
C OUTPUT FILE 'VEL':
C (1) (V(I),I=1,L1234), WHERE L1234 IS THE NUMBER OF GRIDPONTS.
C             L1234=L1*L2*L3*L4.  IF THE FILE 'IND' IS NOT SPECIFIED,
C             L4=N1*N2*N3 BY THE DEFAULT.
C     V(I)... VELOCITY AT THE I-TH GRIDPOINT.
C             FREE SPACE IS INDICATED BY A ZERO VELOCITY V(I)=0.
C
C OUTPUT FILE 'ICB':
C (1) (ICB(I),I=1,L1234), WHERE L1234 IS THE NUMBER OF GRIDPONTS.
C     ICB(I)..INDEX OF (GEOLOGICAL) BLOCK IN WHICH THE I-TH GRIDPOINT
C             IS SITUATED.
C
C OUTPUT FILE 'VELDER':
C (1) FOR EACH GRIDPOINT (1.1):
C (1.1) V,V1,V2,V3,V11,V12,V22,V13,V23,V33
C     VELOCITY AND ITS FIRST AND SECOND DERIVATIVES AT THE GRIDPOINT.
C
C SUBROUTINES REFERENCED:
      EXTERNAL MODEL1,BLOCK,VELOC,PARM2,RARRAY
C     MODEL1,BLOCK,VELOC... FILE 'MODEL.FOR'.
C     PARM2...FILE 'PARM.FOR'.
C     RARRAY..FILE 'ARRAY.FOR'.
C NOTE THAT THE ABOVE SUBROUTINES REFERENCE MANY OTHER EXTERNAL
C PROCEDURES FROM VARIOUS SOURCE CODE FILES OF THE 'MODEL' SUBROUTINE
C PACKAGE.  THESE INDIRECTLY REFERENCED PROCEDURES ARE NOT NAMED HERE,
C BUT ARE LISTED IN THE PARTICULAR SUBROUTINE SOURCE CODE FILES.
C
C DATE: 1994, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     STORAGE LOCATIONS:
C
C     INPUT DATA:
      CHARACTER*80 FGRID,FMODEL,FNET,FIND,FVEL,FICB,FVELD
      INTEGER LU1,LU2,LU3,MPS,MIND
      PARAMETER (LU1=1,LU2=2,LU3=3,MPS=100,MIND=100000)
      INTEGER N1,N2,N3,L1,L2,L3,IPS(MPS),IND(MIND)
      REAL X1MIN,X2MIN,X3MIN,X1MAX,X2MAX,X3MAX
C
C     LU1,LU2,LU3...  LOGICAL UNIT NUMBERS.
C
C     OTHERS:
      LOGICAL LCART,LIND,LVEL,LICB,LVELD,LOBLIQ
      INTEGER MVEL
      PARAMETER (MVEL=10000)
      INTEGER N123,L1234,I1234,IN1,IN2,IN3,IL1,IL2,IL3,IBRICK,INDOLD
      INTEGER ICB(MVEL),NVEL,ISB1,ISRF2,ISB2,ICB2,I,II
      REAL VOUT(MVEL),DX1,DX2,DX3
      REAL COOR(3),UP(10),US(10),VD(10),AUX0,AUX1,AUX2,AUX3,AUX4
      REAL G(12),GAMMA(18),PDER(9),AUX11,AUX12,AUX22,AUX13,AUX23,AUX33
C
C     LCART.. TRUE IF THE MODEL SPECIFIED IN CURVILINEAR COORDINATES IS
C             GRIDDED IN CARTESIAN COORDINATES.
C     LIND... INDICATION OF INDEXED GRID TO SPECIFY IRREGULAR SUBVOLUME
C             OF THE RECTANGULAR VOLUME COVERED BY THE GRID.
C     LVEL,LICB,LVELD... INDICATION OF OPENING AND GENERATING OUTPUT
C             FILES.
C     LOBLIQ..INDICATION OF A SPECIAL OPTION ENABLING TO GRID AN OBLIQUE
C             VERTICAL PROFILE.
C     ICB...  INDICES OF COMPLEX BLOCKS.
C     ISB1... INDEX OF THE SIMPLE BLOCK, SEE SUBROUTINE BLOCK.
C     ISRF2...INDEX OF A SURFACE, SEE SUBROUTINE BLOCK.
C     ISB2... INDEX OF THE SIMPLE BLOCK, SEE SUBROUTINE BLOCK.
C     ICB2... INDEX OF THE COMPLEX BLOCK, SEE SUBROUTINE BLOCK.
C     I...    INDEX OF A GRIDPOINT, OR LOOP VARIABLE.
C     DX1,DX2,DX3... GRID INTERVALS.
C     VEL...  VELOCITY.
C     COOR... CORDINATES OF A GRIDPOINT.
C     UP,US,VD,AUX0,AUX1,AUX2,AUX3,AUX4... AUXILIARY STORAGE LOCATIONS
C             FOR LOCAL MODEL PARAMETERS OR TEMPORARY VARIABLES.
C     G,GAMMA,PDER,AUX11,AUX12,AUX22,AUX13,AUX23,AUX33... AUXILIARY
C             STORAGE LOCATIONS USED IN COORDINATE TRANSFORMATIONS.
C
C.......................................................................
C
      LCART=.TRUE.
C
C.......................................................................
C
C     OPENING FILES AND READING INPUT DATA:
C
C     NAME OF MAIN INPUT DATA:
      WRITE(*,'(A)') ' ENTER THE NAME OF THE MAIN INPUT DATA FILE: '
      FGRID='GRID.DAT'
      READ(*,*) FGRID
      WRITE(*,'(A)') '+                                            '
C
C     READING MAIN INPUT DATA:
      OPEN(LU1,FILE=FGRID,STATUS='OLD')
      FMODEL='MODEL.DAT'
      FNET='NET.DAT'
      FIND=' '
      FVEL='VEL.OUT'
      FICB='ICB.OUT'
      FVELD=' '
      READ(LU1,*) FMODEL,FNET,FIND,FVEL,FICB,FVELD
      DO 11 I=1,MPS
        IPS(I)=0
   11 CONTINUE
      READ(LU1,*) IPS
      DO 15 I=1,MPS
        IF(IPS(I).NE.0) THEN
          IF(IABS(IPS(I)).LT.I) THEN
            IPS(IABS(IPS(I)))=IPS(I)
          ELSE IF(IABS(IPS(I)).GT.I) THEN
            DO 12 II=I+1,MPS
              IF(IPS(II).EQ.0) THEN
                IPS(II)=IPS(I)
                GO TO 13
              END IF
   12       CONTINUE
   13       CONTINUE
          END IF
        END IF
   15 CONTINUE
      DO 16 I=1,MPS
        IF(IPS(I).EQ.0) THEN
          IPS(I)=I
        ELSE IF(IPS(I).GT.0) THEN
          IPS(I)=0
        END IF
   16 CONTINUE
      CLOSE(LU1)
C
C     DATA FOR MODEL:
      OPEN(LU1,FILE=FMODEL,STATUS='OLD')
      CALL MODEL1(LU1)
      CLOSE(LU1)
      IF(KOOR().EQ.0) THEN
C       NO TRANSFORMATION BETWEEN CARTESIAN AND MODEL COORDINATES
        LCART=.FALSE.
      END IF
C
C
C     DATA FOR GRID:
      OPEN(LU1,FILE=FNET,STATUS='OLD')
      N3=1
      L1=1
      L2=1
      L3=1
      READ(LU1,*) N1,N2,N3,L1,L2,L3
      IF(N1.EQ.0) THEN
C       VERTICAL OBLIQUE PROFILE
        IF(L1.NE.0.AND.L1.NE.1) THEN
          PAUSE 'ERROR: INCORRECT L1 FOR AN OBLIQUE PROFILE'
          STOP
        END IF
        LOBLIQ=.TRUE.
        N1=1
        L1=1
      ELSE
        LOBLIQ=.FALSE.
      END IF
      IF(N1.LT.1.OR.N2.LT.1.OR.N3.LT.1.OR.L1.LT.1.OR.L2.LT.1.OR.L3.LT.1)
     *                                                              THEN
        PAUSE 'ERROR: NON-POSITIVE NUMBER OF GRIDPOINTS'
        STOP
      END IF
      X1MIN=0.
      X1MAX=0.
      X2MIN=0.
      X2MAX=0.
      X3MIN=0.
      X3MAX=0.
      READ(LU1,*) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX
      CLOSE(LU1)
C
C     READING THE INDEX ARRAY:
      IF(FIND.EQ.' ') THEN
        LIND=.FALSE.
        IND(1)=1
        L1234=L1*L2*L3*N1*N2*N3
      ELSE
        LIND=.TRUE.
        N123=N1*N2*N3
        IF(N123.GT.MIND) THEN
          PAUSE 'ERROR: TOO MANY BIG BRICKS'
          STOP
        END IF
        DO 31 IBRICK=1,N123
          IND(IBRICK)=IBRICK
   31   CONTINUE
        OPEN(LU1,FILE=FIND,STATUS='OLD')
        READ(LU1,*) (IND(IBRICK),IBRICK=1,N123)
        CLOSE(LU1)
        L1234=0
        DO 32 IBRICK=1,N123
          IF(IND(IBRICK).GT.0) THEN
            L1234=L1234+1
          END IF
   32   CONTINUE
        L1234=L1*L2*L3*L1234
      END IF
C
C     OUTPUT FILE WITH VELOCITIES AT GRIDPOINTS:
      IF(FVEL.EQ.' ') THEN
        LVEL=.FALSE.
      ELSE
        LVEL=.TRUE.
        OPEN(LU1,FILE=FVEL)
      END IF
C
C     OUTPUT FILE WITH INDICES OF COMPLEX BLOCKS:
      IF(FICB.EQ.' ') THEN
        LICB=.FALSE.
      ELSE
        LICB=.TRUE.
        OPEN(LU2,FILE=FICB)
      END IF
C
C     OUTPUT FILE WITH VELOCITY DERIVATIVES AT GRIDPOINTS:
      IF(FVELD.EQ.' ') THEN
        LVELD=.FALSE.
      ELSE
        LVELD=.TRUE.
        OPEN(LU3,FILE=FVELD)
      END IF
C
C.......................................................................
C
C     LOOPS OVER GRIDPOINTS:
C
      IF(LOBLIQ) THEN
        DX1=(X1MAX-X1MIN)/FLOAT(N2*L2)
      ELSE
        DX1=(X1MAX-X1MIN)/FLOAT(N1*L1)
      END IF
      DX2=(X2MAX-X2MIN)/FLOAT(N2*L2)
      DX3=(X3MAX-X3MIN)/FLOAT(N3*L3)
      ISB1=0
      NVEL=0
      IBRICK=0
      INDOLD=0
      I1234=0
      WRITE(*,'(''+'',I18,''  GRIDPOINTS OF'',I8)') I1234,L1234
C
C     LOOP OVER BIG BRICKS:
      DO 83 IN3=0,N3-1
        DO 82 IN2=0,N2-1
          DO 81 IN1=0,N1-1
C
C           CHECK FOR THE COMPUTATIONAL VOLUME:
            IF(LIND) THEN
              IBRICK=IBRICK+1
              IF(IND(IBRICK).EQ.0) THEN
                GO TO 80
              END IF
              IF(IND(IBRICK).NE.INDOLD+1) THEN
                PAUSE 'ERROR: INDICES NOT CONSECUTIVE'
                STOP
              END IF
              INDOLD=IND(IBRICK)
            END IF
C
C           LOOP OVER SMALL BRICKS:
            DO 73 IL3=1,L3
              COOR(3)=X3MIN+DX3*(FLOAT(IN3*L3+IL3)-0.5)
              DO 72 IL2=1,L2
                COOR(2)=X2MIN+DX2*(FLOAT(IN2*L2+IL2)-0.5)
                DO 71 IL1=1,L1
                  IF(LOBLIQ) THEN
                    COOR(1)=X1MIN+DX1*(FLOAT(IN2*L2+IL2)-0.5)
                  ELSE
                    COOR(1)=X1MIN+DX1*(FLOAT(IN1*L1+IL1)-0.5)
                  END IF
C
C                 TRANSFORMATION FROM CARTESIAN TO MODEL COORDINATES:
                  IF(LCART) THEN
                    G(1)=COOR(1)
                    G(2)=COOR(2)
                    G(3)=COOR(3)
                    CALL CARTES(COOR,.FALSE.,G,PDER)
                  END IF
C
C                 VELOCITY EVALUATION:
                  CALL BLOCK(COOR,0,ISB1,ISRF2,ISB2,ICB2,VD)
                  ISB1=ISB2
                  IF(ICB2.EQ.0) THEN
C                   FREE SPACE:
                    DO 41 I=1,10
                      VD(I)=0.
   41               CONTINUE
                  ELSE IF(IPS(ICB2).EQ.0) THEN
C                   DEEMED FREE SPACE:
                    ICB2=0
                    DO 42 I=1,10
                      VD(I)=0.
   42               CONTINUE
                  ELSE
                    CALL PARM2(IABS(ICB2),COOR,UP,US,AUX0,AUX1,AUX2)
                    CALL VELOC(IPS(ICB2),UP,US,
     *                                      AUX1,AUX2,AUX3,AUX4,VD,AUX0)
                  END IF
C
C                 WRITING OUTPUT FILES:
                  IF(NVEL.EQ.MVEL) THEN
                    WRITE(*,'(''+WRITING:  '',I8)') I1234
                    IF(LVEL) THEN
                      CALL WARRAY(LU1,' ','FORMATTED',
     *                                  .FALSE.,0.,.FALSE.,0.,MVEL,VOUT)
C                     FOR VELOCITIES UP TO 9.99999, THE ABOVE STATEMENT
C                     MIGHT BE REPLACED, FOR INSTANCE, BY:
C                     WRITE(LU1,'(10F8.5)') VOUT
                    END IF
                    IF(LICB) THEN
                      WRITE(LU2,'(20(1X,I2))') ICB
                    END IF
                    NVEL=0
                    WRITE(*,'(''+          '')')
                  END IF
                  IF(LVELD) THEN
                    IF(LCART) THEN
C                     TRANSFORMATION FROM MODEL TO CARTESIAN COORDINATES
C                     COVARIANT DERIVATIVES
                      CALL METRIC(COOR,AUX1,G,GAMMA)
                      AUX1=VD(2)
                      AUX2=VD(3)
                      AUX3=VD(4)
                      AUX11=VD( 5)-GAMMA(1)*AUX1-GAMMA( 7)*AUX2
     *                                                   -GAMMA(13)*AUX3
                      AUX12=VD( 6)-GAMMA(2)*AUX1-GAMMA( 8)*AUX2
     *                                                   -GAMMA(14)*AUX3
                      AUX22=VD( 7)-GAMMA(3)*AUX1-GAMMA( 9)*AUX2
     *                                                   -GAMMA(15)*AUX3
                      AUX13=VD( 8)-GAMMA(4)*AUX1-GAMMA(10)*AUX2
     *                                                   -GAMMA(16)*AUX3
                      AUX23=VD( 9)-GAMMA(5)*AUX1-GAMMA(11)*AUX2
     *                                                   -GAMMA(17)*AUX3
                      AUX33=VD(10)-GAMMA(6)*AUX1-GAMMA(12)*AUX2
     *                                                   -GAMMA(18)*AUX3
C                     TRANSFORMATION OF DERIVATIVES
                      VD(2)= AUX1*PDER(1)+ AUX2*PDER(2)+ AUX3*PDER(3)
                      VD(3)= AUX1*PDER(4)+ AUX2*PDER(5)+ AUX3*PDER(6)
                      VD(4)= AUX1*PDER(7)+ AUX2*PDER(8)+ AUX3*PDER(9)
                      AUX1 =AUX11*PDER(1)+AUX12*PDER(2)+AUX13*PDER(3)
                      AUX2 =AUX12*PDER(1)+AUX22*PDER(2)+AUX23*PDER(3)
                      AUX3 =AUX13*PDER(1)+AUX23*PDER(2)+AUX33*PDER(3)
                      VD(5)= AUX1*PDER(1)+ AUX2*PDER(2)+ AUX3*PDER(3)
                      AUX1 =AUX11*PDER(4)+AUX12*PDER(5)+AUX13*PDER(6)
                      AUX2 =AUX12*PDER(4)+AUX22*PDER(5)+AUX23*PDER(6)
                      AUX3 =AUX13*PDER(4)+AUX23*PDER(5)+AUX33*PDER(6)
                      VD(6)= AUX1*PDER(1)+ AUX2*PDER(2)+ AUX3*PDER(3)
                      VD(7)= AUX1*PDER(4)+ AUX2*PDER(5)+ AUX3*PDER(6)
                      AUX1 =AUX11*PDER(7)+AUX12*PDER(8)+AUX13*PDER(9)
                      AUX2 =AUX12*PDER(7)+AUX22*PDER(8)+AUX23*PDER(9)
                      AUX3 =AUX13*PDER(7)+AUX23*PDER(8)+AUX33*PDER(9)
                      VD(8)= AUX1*PDER(1)+ AUX2*PDER(2)+ AUX3*PDER(3)
                      VD(9)= AUX1*PDER(4)+ AUX2*PDER(5)+ AUX3*PDER(6)
                      VD(10)=AUX1*PDER(7)+ AUX2*PDER(8)+ AUX3*PDER(9)
                    END IF
                    WRITE(LU3,'(10F8.5)') VD
                  END IF
                  NVEL=NVEL+1
                  VOUT(NVEL)=VD(1)
                  ICB (NVEL)=ICB2
C
C                 SCREEN OUTPUT:
                  I1234=I1234+1
                  IF(MOD(I1234,1000).EQ.0) THEN
                    WRITE(*,'(''+'',I18)') I1234
                  END IF
C
   71           CONTINUE
   72         CONTINUE
   73       CONTINUE
C
   80       CONTINUE
   81     CONTINUE
   82   CONTINUE
   83 CONTINUE
C
C     REST OF OUTPUT FILES:
      IF(NVEL.GT.0) THEN
        WRITE(*,'(''+WRITING:  '',I8)') I1234
        IF(LVEL) THEN
          CALL WARRAY(LU1,' ','FORMATTED',
     *                                  .FALSE.,0.,.FALSE.,0.,NVEL,VOUT)
C         FOR VELOCITIES UP TO 9.99999, THE ABOVE STATEMENT
C         MIGHT BE REPLACED, FOR INSTANCE, BY:
C         WRITE(LU1,'(10F8.5)') (VOUT(I),I=1,NVEL)
        END IF
        IF(LICB) THEN
          WRITE(LU2,'(20(1X,I2))') (ICB(I),I=1,NVEL)
        END IF
        NVEL=0
      END IF
C
C     SCREEN OUTPUT:
      WRITE(*,'(''+'',I18)') I1234
      STOP
      END
C
C=======================================================================
C
