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