C 
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 Version: 5.30
C Date: 1999, March 28
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
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 Subroutine CARTES
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.......................................................................
C
C                                                    
C Description of data files:
C
C Name of the main input data file read from the * external unit:
C (1) 'SEP','GRID'
C     'SEP'...String in apostrophes containing the name of the input
C             file with the data specifying grid dimensions and the
C             name of the data specifying the model.
C             Description of file SEP
C     'GRID'..Name of the input data file described below.
C     Defaults: 'SEP'='grid.h', 'GRID'='grid.dat'.
C
C Main input data file GRID:
C (1) 'IND','VEL','ICB','VEL1','VEL2','VEL3','VEL11',
C     'VEL12','VEL22','VEL13','VEL23','VEL33',/
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     'VEL1','VEL2','VEL3','VEL11','VEL12','VEL22','VEL13','VEL23',
C     'VEL33'... Names of the output formatted files containing
C             individual first or second partial velocity derivatives
C             at the given gridpoints.
C             If the filename is blank, the corresponding file is not
C             created.
C     Defaults: 'IND'=' ', 'VEL'='vel.out', 'ICB'='icb.out', 'VEL1'=' ',
C             'VEL2'=' ', 'VEL3'=' ', 'VEL11'=' ', 'VEL12'=' ',
C             'VEL22'=' ', 'VEL13'=' ', 'VEL23'=' ', 'VEL33'=' '.
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 Example of data set GRID
C
C                                                     
C Data file SEP has the form of the SEP (Stanford Exploration Project)
C parameter file:
C     All the data are specified in the form of PARAMETER=VALUE, e.g.
C     N1=50, with PARAMETER directly preceding = without intervening
C     spaces and with VALUE directly following = without intervening
C     spaces.  The PARAMETER=VALUE couple must be delimited by a space
C     or comma from both sides.
C     The PARAMETER string is not case-sensitive.
C     PARAMETER= followed by a space resets the default parameter value.
C     All other text in the input files is ignored.  The file thus may
C     contain unused data or comments without leading comment character.
C     Everything between comment character # and the end of the
C     respective line is ignored, too.
C     The PARAMETER=VALUE couples may be specified in any order.
C     The last appearance takes precedence.
C Input file specifying the model:
C     MODEL='string'... Input data file describing the model, it is
C             described in the subroutine file 'model.for'.
C             Description of file MODEL
C             Default:  'MODEL'='model.dat'
C Data specifying grid dimensions:
C     N1=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1=1
C             Special option of N1=0:
C               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               Here
C               X1MIN=O1-0.5*D1,  X1MAX=X1MIN+FLOAT(N1)*D1,
C               X2MIN=O2-0.5*D2,  X2MAX=X2MIN+FLOAT(N2)*D2,
C               X3MIN=O3-0.5*D3,  X3MAX=X3MIN+FLOAT(N3)*D3.
C     N2=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2=1
C     N3=positive integer... Number of gridpoints along the X3 axis.
C             Default: N3=1
C     D1=real... Grid interval in the direction of the first coordinate
C             axis.
C             Default: D1=1.
C     D2=real... Grid interval in the direction of the second coordinate
C             axis.
C             Default: D2=1.
C     D3=real... Grid interval in the direction of the third coordinate
C             axis.
C             Default: D3=1.
C     O1=real... First coordinate of the grid origin (first point of the
C             grid).
C             Default: O1=0.
C     O2=real... Second coordinate of the grid origin.
C             Default: O2=0.
C     O3=real... Third coordinate of the grid origin.
C             Default: O3=0.
C Additional parameters for 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.  Here
C               X1MIN=O1-0.5*D1,  X1MAX=X1MIN+FLOAT(N1)*D1,
C               X2MIN=O2-0.5*D2,  X2MAX=X2MIN+FLOAT(N2)*D2,
C               X3MIN=O3-0.5*D3,  X3MAX=X3MIN+FLOAT(N3)*D3.
C             Then (if the numbers L1,L2,L3 are specified in addition to
C             N1,N2,N3) each big brick is subdivided into L1*L2*L3 small
C             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             A special option of N1=0:
C               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
C               cells, each big cell being divided into L2*L3 small
C               cells.
C     L1=positive integer... Number of small bricks in one big brick in
C             the direction of axis X1.  If specified, must be positive.
C             Default: L1=1
C     L2=positive integer... Number of small bricks in one big brick in
C             the direction of axis X2.  If specified, must be positive.
C             Default: L2=1
C     L3=positive integer... Number of small bricks in one big brick in
C             the direction of axis X3.  If specified, must be positive.
C             Default: L3=1
C             (If the numbers L1,L2,L3 are not specified, each big brick
C             contains just one small brick, as large as big one.)
C Example of data set SEP
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 files 'VEL','VEL1','VEL2','VEL3','VEL11','VEL12','VEL22',
C 'VEL13','VEL23':
C (1) (V(I),I=1,L1234), where L1234 is the number of gridpoints.
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 or its partial derivative at the I-th gridpoint.
C             Free space is indicated by a zero velocity or derivative
C             V(I)=0.
C
C Output file 'ICB':
C (1) (ICB(I),I=1,L1234), where L1234 is the number of gridpoints.
C     ICB(I)..Index of (geological) block in which the I-th gridpoint
C             is situated.
C
C.......................................................................
C
C Subroutines referenced:
      EXTERNAL KOOR,MODEL1,BLOCK,VELOC,PARM2,WARRAY
      INTEGER  KOOR
C     KOOR... File 'metric.for'.
C     MODEL1,BLOCK,VELOC... File 'model.for'.
C     PARM2...File 'parm.for'.
C     WARRAY..File 'forms.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-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C     Allocation of arrays:
      INTEGER MVEL,MIND
      PARAMETER (MVEL=10000,MIND=MRAM-11*MVEL)
      INTEGER ICB(MVEL),IND(MIND)
      REAL VOUT(MVEL),VOUD(MVEL,9)
      EQUIVALENCE (VOUT,RAM)
      EQUIVALENCE (ICB,RAM(MVEL+1))
      EQUIVALENCE (VOUD,RAM(2*MVEL+1))
      EQUIVALENCE (IND,RAM(11*MVEL+1))
C
C.......................................................................
C
C     Storage locations:
C
C     Input data:
      CHARACTER*80 FGRID,FMODEL,FSEP,FIND,FVEL,FICB,FVELD(9)
      INTEGER LU1,LU2,LUD(9),MPS
      PARAMETER (LU1=1,LU2=2,MPS=100)
      INTEGER N1,N2,N3,L1,L2,L3,IPS(MPS)
      REAL D1,D2,D3,O1,O2,O3
      REAL X1MIN,X2MIN,X3MIN,X1MAX,X2MAX,X3MAX
C
C     LU1,LU2,LUD...  Logical unit numbers.
C
C     Others:
      LOGICAL LCART,LIND,LVEL0,LICB,LVELD,LVEL(9),LOBLIQ
      INTEGER N123,L1234,I1234,IN1,IN2,IN3,IL1,IL2,IL3,IBRICK,INDOLD
      INTEGER NVEL,ISRF2,ISB2,ICB2,I,II
      REAL 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     LVEL0,LICB,LVELD,LVEL... Indication of opening and generating
C             output files.
C     LOBLIQ..Indication of a special option enabling to grid an oblique
C             vertical profile.
C     ICB...  Indices of complex blocks.
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... Coordinates 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
      DATA LUD/3,4,5,6,7,8,9,10,11/
C
C.......................................................................
C
      LCART=.TRUE.
C
C.......................................................................
C
C     Opening files and reading input data:
C
C     Name of main input data:
      FSEP='grid.h'
      FGRID='grid.dat'
      WRITE(*,'(A)') '+GRID: Enter input filenames [grid.h, grid.dat]: '
      READ(*,*) FSEP,FGRID
      WRITE(*,'(A)') '+GRID: Reading data...                           '
      CALL RSEP1(LU1,FSEP)
C
C     Reading main input data:
      OPEN(LU1,FILE=FGRID,STATUS='OLD')
      FIND=' '
      FVEL='vel.out'
      FICB='icb.out'
      DO 10 I=1,9
        FVELD(I)=' '
   10 CONTINUE
      READ(LU1,*) FIND,FVEL,FICB,(FVELD(I),I=1,9)
      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=IABS(IPS(I)),MPS
              IF(IPS(II).EQ.0) THEN
                IPS(II)=IPS(I)
                IPS(I)=0
                GO TO 13
              END IF
   12       CONTINUE
C           GRID-01
            CALL ERROR('GRID-01: Too small array IPS(MPS)')
   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:
      CALL RSEP3T('MODEL',FMODEL,'model.dat')
      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     Data for grid:
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3I('L1',L1,1)
      CALL RSEP3I('L2',L2,1)
      CALL RSEP3I('L3',L3,1)
      IF(N1.EQ.0) THEN
C       Vertical oblique profile
        IF(L1.NE.0.AND.L1.NE.1) THEN
C         GRID-02
          CALL ERROR('GRID-02: Incorrect L1 for an oblique profile')
        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
C       GRID-03
        CALL ERROR('GRID-03: Non-positive number of gridpoints')
      END IF
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
      CALL RSEP3R('O1',O1,0.)
      CALL RSEP3R('O2',O2,0.)
      CALL RSEP3R('O3',O3,0.)
      X1MIN=O1-0.5*D1
      X2MIN=O2-0.5*D2
      X3MIN=O3-0.5*D3
      X1MAX=X1MIN+FLOAT(N1)*D1
      X2MAX=X2MIN+FLOAT(N2)*D2
      X3MAX=X3MIN+FLOAT(N3)*D3
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
C         GRID-04
          CALL ERROR('GRID-04: Too many big bricks')
        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
        LVEL0=.FALSE.
      ELSE
        LVEL0=.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:
      LVELD=.FALSE.
      DO 33 I=1,9
        IF(FVELD(I).EQ.' ') THEN
          LVEL(I)=.FALSE.
        ELSE
          LVELD  =.TRUE.
          LVEL(I)=.TRUE.
          OPEN(LUD(I),FILE=FVELD(I))
        END IF
   33 CONTINUE
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)
      NVEL=0
      IBRICK=0
      INDOLD=0
      I1234=0
      WRITE(*,'(''+GRID: '',I16,''  gridpoints of'',I9)') 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
C               GRID-05
                CALL ERROR('GRID-05: Indices not consecutive')
              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,0,ISRF2,ISB2,ICB2)
                  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(*,'(''+GRID: Writing'',I9)') I1234
                    IF(LVEL0) 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
                    DO 51 I=1,9
                      IF(LVEL(I)) THEN
                        CALL WARRAY(LUD(I),' ','FORMATTED',
     *                             .FALSE.,0.,.FALSE.,0.,MVEL,VOUD(1,I))
                      END IF
   51               CONTINUE
                    NVEL=0
                    WRITE(*,'(''+GRID:        '')')
                  END IF
                  NVEL=NVEL+1
                  VOUT(NVEL)=VD(1)
                  ICB (NVEL)=ICB2
                  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
                    DO 61 I=1,9
                      IF(LVEL(I)) THEN
                        VOUD(NVEL,I)=VD(I+1)
                      END IF
   61               CONTINUE
                  END IF
C
C                 Screen output:
                  I1234=I1234+1
                  IF(MOD(I1234,1000).EQ.0) THEN
                    WRITE(*,'(''+GRID: '',I16)') 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(*,'(''+GRID: Writing'',I9)') I1234
        IF(LVEL0) 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
        DO 91 I=1,9
          IF(LVEL(I)) THEN
            CALL WARRAY(LUD(I),' ','FORMATTED',
     *                             .FALSE.,0.,.FALSE.,0.,NVEL,VOUD(1,I))
          END IF
   91   CONTINUE
        NVEL=0
      END IF
C
C     Screen output:
      WRITE(*,'(''+GRID: done'',I12)') I1234
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'model.for'
C     model.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'srfc.for'
C     srfc.for
      INCLUDE 'parm.for'
C     parm.for
      INCLUDE 'val.for'
C     val.for
      INCLUDE 'fit.for'
C     fit.for
C
C=======================================================================
C