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