C 
C Program GRID to discretize functions, specifying velocities and other
C material parameters or describing structural interfaces, at gridpoints
C of a regular rectangular grid.
C
C Useful for the full wave finite differences, the shortest path
C calculation of seismic rays, eikonal equation 'finite differences',
C raster imaging of the model, or generating test data for inversion.
C Note that also an oblique vertical 2-D section across the 3-D model
C may be gridded.
C
C Version: 5.40
C Date: 2000, May 15
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 value of the
C input parameter KOORGRD=0 (default) to KOORGRD=1.
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
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     KOORGRD=integer... Coordinates for discretization:
C             KOORGRD=0: Cartesian coordinates.
C             KOORGRD=1: Model coordinates.
C             Default: KOORGRD=0
C Data specifying the function to be gridded:
C     ISRF=integer...
C             ISRF=0: Material parameters in complex blocks will be
C               gridded (default).
C             ISRF=positive integer:  Function describing surface ISRF
C               will be gridded.
C             Default: ISRF=0
C     ICBEXT=integer... Used if ISRF=0:
C             ICBEXT=0: Indices of complex blocks will be determined at
C               the gridpoints and the material parameters will
C               correspond to the complex blocks (default).
C             ICBEXT=positive integer:  Material parameters of complex
C               block ICBEXT will be calculated at all gridpoints
C               (complex block ICBEXT is extended to the whole grid).
C             Default: ICBEXT=0
C     MPAR=integer... Material parameter to be gridded.
C             Used if ISRF=0 and ICBEXT=0:
C             MPAR=0: The material parameter for each block is listed in
C               input file LPAR (default).
C             MPAR=positive integer: Parameter number MPAR will be
C               gridded in each complex block:
C               MPAR=1: P wave velocity,
C               MPAR=2: S wave velocity,
C               MPAR=3: density,
C               MPAR=4: P wave loss factor,
C               MPAR=5: S wave loss factor,
C               MPAR=6 to 26: Reduced (i.e., divided by the density)
C                 anisotropic elastic parameters A11, A12, A22, A13,
C                 A23, A33, A14, A24, A34, A44, A15, A25, A35, A45, A55,
C                 A16, A26, A36, A46, A56 or A66.
C               MPAR=27-47: Reduced (i.e., divided by the density)
C                 imaginary parts of anisotropic elastic parameters Q11,
C                 Q12, Q22, Q13, Q23, Q33, Q14, Q24, Q34, Q44, Q15, Q25,
C                 Q35, Q45, Q55, Q16, Q26, Q36, Q46, Q56 or Q66.
C             Default: MPAR=0
C     LPAR='string'... Name of the input formatted file containing the
C             list of material parameters to be gridded.
C             The file should contain one integer MPAR(k) per each
C             complex block k.  The meaning of integers MPAR(k) is
C             analogous to input parameter MPAR.  In addition, MPAR(k)=0
C             means zeros in the corresponding complex block.
C             The default values of all integers MPAR(k) are 1 (P-wave
C             velocity).  LPAR=' ' (default) has the same meaning as
C             MPAR=1.
C             Used if ISRF=0 and ICBEXT=0 and MPAR=0.
C             Default: LPAR=' '
C Data specifying filenames with gridded values:
C     IND='string'... 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             Description of file IND
C             Default: IND=' '
C     ICB='string'... Name of the output formatted file containing the
C             indices of complex geological blocks at the gridpoints.
C             If blank, the file is not created.
C             Description of the output files
C             Default: ICB=' '
C     VEL='string'... Name of the output formatted file,
C             containing the velocities, other material parameters or
C             the values of the function describing the given surface
C             at the gridpoints.
C             Velocity=0 is inserted in a free space.
C             If blank, the file is not created.
C             Description of the output files
C             Default: VEL='vel.out'
C     VEL1='string', VEL2='string', VEL3='string', VEL11='string',
C     VEL12='string', VEL22='string', VEL13='string', VEL23='string',
C     VEL33='string'... 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             Description of the output files
C             Defaults: VEL1=' ', VEL2=' ', VEL3=' ', VEL11=' ',
C             VEL12=' ', VEL22=' ', VEL13=' ', VEL23=' ', VEL33=' '
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                                                     
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                                                     
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 For general description of the files with gridded data refer to file
C forms.htm.
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 /MODELC/:
      INCLUDE 'model.inc'
C     model.inc
C None of the storage locations of the common block are altered.
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 ISRF,ICBEXT,MPAR
      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 LIND,LVEL0,LICB,LVELD,LVEL(9),LOBLIQ
      INTEGER KOORG
      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
      REAL A(10,21),Q(21)
C
C     KOORG...0 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
C     Opening files and reading input data:
C
C     Name of main input data:
      FSEP=' '
      WRITE(*,'(A)') '+GRID: Enter input filename: '
      READ(*,*) FSEP
C
C     Reading all data from the SEP file into the memory:
      IF (FSEP.NE.' ') THEN
        CALL RSEP1(LU1,FSEP)
      ELSE
C       GRID-07
        CALL ERROR('GRID-07: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
      WRITE(*,'(A)') '+GRID: Reading data...       '
C
C     Reading filenames:
      CALL RSEP3T('IND'  ,FIND    ,' ')
      CALL RSEP3T('ICB'  ,FICB    ,' ')
      CALL RSEP3T('VEL'  ,FVEL    ,'vel.out')
      CALL RSEP3T('VEL1' ,FVELD(1),' ')
      CALL RSEP3T('VEL2' ,FVELD(2),' ')
      CALL RSEP3T('VEL3' ,FVELD(3),' ')
      CALL RSEP3T('VEL11',FVELD(4),' ')
      CALL RSEP3T('VEL12',FVELD(5),' ')
      CALL RSEP3T('VEL22',FVELD(6),' ')
      CALL RSEP3T('VEL13',FVELD(7),' ')
      CALL RSEP3T('VEL23',FVELD(8),' ')
      CALL RSEP3T('VEL33',FVELD(9),' ')
C
C     Data for model:
      CALL RSEP3T('MODEL',FMODEL,'model.dat')
      OPEN(LU1,FILE=FMODEL,STATUS='OLD')
      CALL MODEL1(LU1)
      CLOSE(LU1)
      CALL RSEP3I('KOORGRD',KOORG,0)
      IF(KOOR().EQ.0) THEN
C       No transformation between Cartesian and model coordinates
        KOORG=1
      END IF
C
C     Reading indices of material parameters:
      CALL RSEP3I('ISRF',ISRF,0)
      CALL RSEP3I('ICBEXT',ICBEXT,0)
      CALL RSEP3I('MPAR',MPAR,0)
      IF(MPS.LT.NCB) THEN
C       GRID-01
        CALL ERROR('GRID-01: Too small array IPS(MPS)')
      END IF
      DO 11 I=1,NCB
        IPS(I)=MAX0(1,MPAR)
   11 CONTINUE
      IF(MPAR.EQ.0) THEN
        CALL RSEP3T('LPAR',FGRID,' ')
        IF(FGRID.NE.' ') THEN
          OPEN(LU1,FILE=FGRID,STATUS='OLD')
          READ(LU1,*) (IPS(I),I=1,NCB)
          CLOSE(LU1)
        END IF
      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(KOORG.EQ.0) 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:
                  IF(ISRF.EQ.0) THEN
                    IF(ICBEXT.EQ.0) THEN
                      CALL BLOCK(COOR,0,0,ISRF2,ISB2,ICB2)
                    ELSE
                      ICB2=ICBEXT
                    END IF
                    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 IF(IPS(ICB2).LE.5) THEN
C                     Isotropic elastic parameters:
                      CALL PARM2(IABS(ICB2),COOR,UP,US,AUX0,AUX1,AUX2)
                      IF(IPS(ICB2).EQ.1.OR.IPS(ICB2).EQ.4) THEN
                        CALL VELOC( 1,UP,US,AUX1,AUX2,AUX3,AUX4,VD,AUX0)
                      ELSE IF(IPS(ICB2).EQ.2.OR.IPS(ICB2).EQ.5) THEN
                        CALL VELOC(-1,UP,US,AUX1,AUX2,AUX3,AUX4,VD,AUX0)
                      END IF
                      IF(IPS(ICB2).GT.2) THEN
                        VD(1)=AUX0
                        DO 43 I=2,10
                          VD(I)=0.
   43                   CONTINUE
                      END IF
                    ELSE IF(IPS(ICB2).LE.47) THEN
C                     Anisotropic elastic parameters:
                      CALL PARM3(IABS(ICB2),COOR,A,AUX0,Q)
                      IF(IPS(ICB2).LE.26) THEN
                        DO 46 I=1,10
                          VD(I)=A(I,IPS(ICB2)-5)
   46                   CONTINUE
                      ELSE
                        VD(1)=Q(IPS(ICB2)-26)
                        DO 47 I=2,10
                          VD(I)=0.
   47                   CONTINUE
                      END IF
                    ELSE
C                     
C                     GRID-06
                      CALL ERROR('GRID-06: Wrong material parameter')
C                     Material parameters are indexed by integers
C                     from 1 to 47, but index greater than 47 has been
C                     encountered.  Check the input data.
                    END IF
                  ELSE
                    CALL SRFC2(ISRF,COOR,VD)
                  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(KOORG.EQ.0) 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