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: 6.60 C Date: 2012, February 10 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/klimes.htm 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 C 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 NEGPAR=integer... Flag whether the negative values of material C parameters are allowed. For the description refer to file C model.for. 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: 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 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 if C ICBEXT=0. Otherwise, the file would be filled by the C value of ICBEXT. C If ICB is blank (default), 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(N2)*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 Optional parameters specifying the form of the quantities C written in the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C C forms.for. C NUMLIN=positive integer ... Number of reals or integers in one C line of the output file. See the description in file C C forms.for. C Optional parameter for screen output control: C NWRITE=positive integer ... Screen output is updated after C NWRITE gridpoints are calculated. C Default: NWRITE=1000 C Example of data set C 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),NWRITE 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.) CALL RSEP3I('NWRITE',NWRITE,1000) X1MIN=O1-0.5*D1 X2MIN=O2-0.5*D2 X3MIN=O3-0.5*D3 IF(LOBLIQ) THEN X1MAX=X1MIN+FLOAT(N2)*D1 ELSE X1MAX=X1MIN+FLOAT(N1)*D1 END IF 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 C 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 CALL WARRAI(LU2,' ','FORMATTED', * .FALSE.,0,.FALSE.,0,MVEL,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,NWRITE).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 CALL WARRAI(LU2,' ','FORMATTED',.FALSE.,0,.FALSE.,0,NVEL,ICB) 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