C
C Provisional version MODLE2D of the program designed to calculate C directional Lyapunov exponents and the mean Lyapunov exponent C for a 2-D model without interfaces. 3-D models may be used if NY C is specified. C C Version: 5.70 C Date: 2003, May 12 C C Coded by Ludek Klimes C klimes@seis.karlov.mff.cuni.cz 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 Data specifying the model: C MODEL='string'... String containing the name of the input data C file specifying the model. C Description of file MODEL C No default, MODEL must be specified and cannot be blank. C Data to specify the calculation of Lyapunov exponents: C KOOR1=integer... Index of the first coordinate determining the C 2-D section for the calculation of the Lyapunov exponents. C KOOR1 must be 1, 2 or 3. C Default: KOOR1=1 C KOOR2=integer... Index of the second coordinate determining the C 2-D section for the calculation of the Lyapunov exponents. C KOOR2 must be 1, 2 or 3 and must differ from KOOR1. C Default: KOOR2=2 C NA=integer... Number of angular directions for the calculation of C the directional Lyapunov exponents. C Default: NA=90 C DA... Angular step for the calculation of the directional C Lyapunov exponents. The default value is recommended. C Default: DA=3.141592/NA C OA... The first angle for the calculation of the directional C Lyapunov exponents. The angles are defined in radians, C -pi/2 for positive half-axis KOOR1, 0 for positive C half-axis KOOR2, pi/2 for negative half-axis KOOR1. C The default value is usually sufficient. C Default: OA=-1.570796+0.5*DA C NY=integer... Number of sections for the calculation of each C directional Lyapunov exponent. The sections are C regularly spaced through the model volume. C Default: NY=1 C NX=integer... Number of straight lines for the calculation of C each directional Lyapunov exponent. The lines are C equally spaced and cover the 2-D section of the model C box. C Default: NX=45 C NS=integer... Number of segments along the longest line for the C calculation of a directional Lyapunov exponent. The C length of the segments determined from NS is then used C for all lines corresponding to the direction. C Default: NS=45 C Name of the optional input file: C MODLEW='string' ... Name of the optional input file with weights C for individual angular directions. If specified, the C directional Lyapunov exponents are multiplied by C the weights. The weights thus affect also the mean C Lyapunov exponent for the model. The maximum directional C Lyapunov exponent of the model is not affected. C Description of file MODLEW C Default: MODLEW=' ' (No weights) C Data to specify the frame of the graph of the Lyapunov exponents: C LEMAX...Value of the Lyapunov exponent corresponding to the top C of the frame. The bottom of the frame cooresponds to 0. C Default: LEMAX=(maximum directional Lyapunov exponent) C ALEMIN..Angle corresponding to the left-hand edge of the frame. C The default value is usually sufficient. C Default: ALEMIN=OA-0.5*DA C ALEMAX..Angle corresponding to the right-hand edge of the frame. C The default value is usually sufficient. C Default: ALEMAX=OA+DA*(NA-0.5) C Names of the output files: C MODLED='string'... Name of the output file containing, for each C angle, the directional Lyapunov exponent, optionally C multiplied by the corresponding weight, if the input file C MODLEW is given. The file has form C LINes to enable C plotting by program C pictures.for. C The first coordinate of a point of a line is the angle, C the second coordinate is the corresponding directional C Lyapunov exponent. The name of the line is the value of C the maximum directional Lyapunov exponent and the C reference point of the line should enable to draw this C value. C Default: MODLED='modled.out' C MODLEM='string'... Name of the output file containing the mean C Lyapunov exponent for the model (a single value). C The mean Lyapunov exponent is calculated by averaging C directional Lyapunov exponents. The file has form C LINes to enable C to plot the mean Lyapunov exponent as a horizontal line C into the graph of the directional Lyapunov exponents. C The line has two points, the first coordinate is the C minimum or maximum angle, respectively. The second C coordinate is the mean Lyapunov exponent. The name of C the line is the value of the mean Lyapunov exponent and C the reference point of the line should enable to draw C this value. C Default: MODLEM='modlem.out' C MODLEF='string'... Name of the output file containing the lines C composing the frame of the graph of the directional C Lyapunov exponents. The file has form C LINes, similarly C as files MODLED and MODLEM. C Default: MODLEF='modlef.out' C C C Optional input data file MODLEW: C This data file contains the weights of directional Lyapunov C exponents. All NA weights must be specified. Each record 2.2 C (see below) containes the angle A and the corresponding weight W. C The angles must be in agreement with the input parameters OA, DA. C The file has form C LINes. C The data are read in by the list directed input (free format). C In the list of input data below, each numbered paragraph C indicates the beginning of a new input operation (new READ C statement). C (1) None to several (20 at most) strings terminated by / (a slash). C (2) The weights - the header 2.1 and NA times the values 2.2: C (2.1) 'Text',A,W,/ ... Reference text and optional reference C point specified as angle and weight. Not used in present C version, / (a slash) only is thus sufficient. C (2.2) A,W,/ ... Angle A and the corresponding weight W. The data 2.2 C must repeat NA times. C A ... Angle of the direction, must correspond to the input C parameters OA,DA, i.e. I-th angle must equal OA+(I-1)*DA. C W ... Weight corresponding to the direction A. C (3) / (a slash) or the end of file. C----------------------------------------------------------------------- C C Common block /MODELC/: INCLUDE 'model.inc' C None of the storage locations of the common block are altered. C C----------------------------------------------------------------------- C CHARACTER*80 FILE1,FILE2 CHARACTER*1 CH(20) CHARACTER*30 TEXT INTEGER LU1,LU2 PARAMETER (LU1=1,LU2=2) LOGICAL LINEND INTEGER KOOR1,KOOR2,KOOR3,NA,NX,NY,NS,IA,IX,IY,IS,NUMSUM,KS,I, * ISRF,ISB,ICB REAL COOR(3),H(3) REAL UP(10),US(10),VD(10),AUX,AUX0,AUX1,AUX2,AUX3,AUX4,OA,DA,A,O1, * O2,D1,D2,DX,DS0,DY,BOUND1,BOUND2,BOUND3,BOUND4,BOUND5,BOUND6, * EL2MOD,EL2MAX,EL2SUM,TTSUM,X1,X2,S,DS,TT,ZON,ELA,EL0,EL1,EL2, * ELAOLD,PHI,PHIOLD,PHI0,COOR1,COOR2,VV,V1,V2,V3,X1OLD,X2OLD, * VOLD,VVOLD,VVNEG,VVPOS,AW,W,EMAX,ALEMIN,ALEMAX C C....................................................................... C C Main input data: WRITE(*,'(A)') '+MODLE2D: Enter input filename: ' FILE1=' ' READ(*,*) FILE1 IF(FILE1.EQ.' ') THEN C MODLE2D-01 CALL ERROR('MODLE2D-01: No input SEP file specified') END IF CALL RSEP1(LU1,FILE1) WRITE(*,'(A)') '+MODLE2D: Working... ' C C Data for model: CALL RSEP3T('MODEL',FILE1,' ') IF(FILE1.EQ.' ') THEN C MODLE2D-02 CALL ERROR('MODLE2D-02: No model specified') END IF OPEN(LU1,FILE=FILE1,STATUS='OLD') CALL MODEL1(LU1) CLOSE(LU1) C C Input parameters CALL RSEP3I('KOOR1',KOOR1,1) CALL RSEP3I('KOOR2',KOOR2,2) CALL RSEP3I('NA',NA,90) CALL RSEP3R('DA',DA,3.141592/FLOAT(NA)) CALL RSEP3R('OA',OA,-1.570796+0.5*DA) IF(OA.LT.-1.570796.OR.OA+DA*FLOAT(NA-1).GT.1.570796) THEN C MODLE2D-03 CALL ERROR('MODLE2D-03: Wrong direction') END IF C Number of sections CALL RSEP3I('NY',NY,1) C Number of lines CALL RSEP3I('NX',NX,45) C Maximum number of steps along a line CALL RSEP3I('NS',NS,45) C C 2-D section BOUND1=BOUNDM(2*KOOR1-1) BOUND2=BOUNDM(2*KOOR1) BOUND3=BOUNDM(2*KOOR2-1) BOUND4=BOUNDM(2*KOOR2) KOOR3=6-KOOR2-KOOR1 BOUND5=BOUNDM(2*KOOR3-1) BOUND6=BOUNDM(2*KOOR3) C C Output file for directional Lyapunov exponents: CALL RSEP3T('MODLED',FILE1,'modled.out') IF(FILE1.NE.' ') THEN OPEN(LU1,FILE=FILE1) WRITE(LU1,'(A)') '/' END IF C C Optional file with the weights: CALL RSEP3T('MODLEW',FILE2,' ') IF (FILE2.NE.' ') THEN OPEN(LU2,FILE=FILE2) READ(LU2,*) CH READ(LU2,*) CH(1) ENDIF C C Loop over directions EL2MOD=0. EL2MAX=0. WRITE(*,'(3(A,I3))') '+MODLE2D:',0,' /',NA,' directions' DO 90 IA=0,NA-1 A=OA+DA*FLOAT(IA) C A is the deviation from the X2 axis in radians O1=BOUND1 O2=BOUND3 D1=BOUND2-BOUND1 D2=BOUND4-BOUND3 DX =(ABS(D1*COS(A))+ABS(D2*SIN(A)))/FLOAT(NX) DS0=(ABS(D1*SIN(A))+ABS(D2*COS(A)))/FLOAT(NS) C *** Initiating averaging over lines NUMSUM=0 TTSUM =0. * ZONSUM=0. * ELASUM=0. * PHISUM=0. * EL1SUM=0. EL2SUM=0. C *** End of initialization DY=(BOUND6-BOUND5)/FLOAT(NY) C Loop over sections DO 85, IY=0,NY-1 COOR(KOOR3)=BOUND5+DY/2.+FLOAT(IY)*DY C Loop over lines DO 80 IX=0,NX-1 C Initial point of the line KS=0 IF(SIN(A).LE.0.) THEN X1=O1+(FLOAT(IX)+0.5)*DX/COS(A) X2=O2 IF(X1.GT.BOUND2) THEN X2=O2+(BOUND2-X1)*COS(A)/SIN(A) X1=BOUND2 END IF ELSE X1=BOUND2-(FLOAT(IX)+0.5)*DX/COS(A) X2=O2 IF(X1.LT.BOUND1) THEN X2=O2+(BOUND1-X1)*COS(A)/SIN(A) X1=BOUND1 END IF END IF C Unit vetor perpendicular to the line H(KOOR1)=COS(A) H(KOOR2)=-SIN(A) H(6-KOOR1-KOOR2)=0. C C *** Initiating numerical quadrature along a part of the line 20 CONTINUE S =0. TT =0. ZON=0. ELA=0. PHI=0. EL0=0. EL1=0. EL2=0. ELAOLD=0. PHIOLD=0. PHI0=0. C *** End of initialization LINEND=.FALSE. C Loop over points of a line DO 70 IS=KS,NS+1 COOR1=X1+DS0*SIN(A)*FLOAT(IS) COOR2=X2+DS0*COS(A)*FLOAT(IS) IF(COOR1.LT.BOUND1) THEN COOR1=BOUND1 COOR2=X2+(COOR1-X1)*COS(A)/SIN(A) LINEND=.TRUE. END IF IF(COOR1.GT.BOUND2) THEN COOR1=BOUND2 COOR2=X2+(COOR1-X1)*COS(A)/SIN(A) LINEND=.TRUE. END IF IF(COOR2.GE.BOUND4) THEN COOR2=BOUND4 COOR1=X1+(COOR2-X2)*SIN(A)/COS(A) LINEND=.TRUE. END IF C COOR1,COOR2 is the current point on the line C *** Beginning of numerical quadrature COOR(KOOR1)=COOR1 COOR(KOOR2)=COOR2 CALL BLOCK(COOR,0,0,ISRF,ISB,ICB) IF(ICB.EQ.0) THEN C Free space VD(1)=0. VV=0. ELSE C Material complex block CALL PARM2(IABS(ICB),COOR,UP,US,AUX0,AUX1,AUX2) CALL VELOC(ICB,UP,US,AUX1,AUX2,AUX3,AUX4,VD,AUX0) V1=VD( 5)*H(1)+VD( 6)*H(2)+VD( 8)*H(3) V2=VD( 6)*H(1)+VD( 7)*H(2)+VD( 9)*H(3) V3=VD( 8)*H(1)+VD( 9)*H(2)+VD(10)*H(3) VV=V1 *H(1)+V2 *H(2)+V3 *H(3) VV=VV/VD(1) IF(IS.GT.KS) THEN DS=SQRT((COOR1-X1OLD)**2+(COOR2-X2OLD)**2) S=S+DS TT=TT+DS*(0.5/VD(1)+0.5/VOLD) IF(VVOLD*VV.GE.0.) THEN ELA=ELA+DS*VVNEG/2. PHI=PHI+DS*VVPOS/2. ELSE ELA=ELA+DS*VVNEG*ABS(VVOLD/(VVOLD-VV))*2./3. PHI=PHI+DS*VVPOS*ABS(VVOLD/(VVOLD-VV))*2./3. END IF IF(VVOLD.LT.0..AND.VV.GE.0.) THEN C End of negative VV ("high velocity") IF(EL0.EQ.0.) THEN C First end of negative VV ("high velocity") EL1=EL1+AMAX1(0.,ELA-ELAOLD) EL2=EL2+AMAX1(0.,ELA-ELAOLD) PHI0=AMIN1(PHI-PHIOLD,0.523599) ELSE EL1=EL1 * +AMAX1(0.,ELA-ELAOLD+ALOG(ABS(COS(PHI-PHIOLD)))) EL2=EL2+AMAX1(0.,ELA-ELAOLD) * +ALOG(COS(AMIN1(PHI-PHIOLD,1.047198))) END IF ZON=ZON+1. EL0=EL0+AMAX1(0.,ELA-ELAOLD) ELAOLD=ELA PHIOLD=PHI END IF VVNEG=SQRT(AMAX1(0.,-VV)) VVPOS=SQRT(AMAX1(0., VV)) IF(VVOLD*VV.GE.0.) THEN ELA=ELA+DS*VVNEG/2. PHI=PHI+DS*VVPOS/2. ELSE ELA=ELA+DS*VVNEG*ABS(VV/(VVOLD-VV))*2./3. PHI=PHI+DS*VVPOS*ABS(VV/(VVOLD-VV))*2./3. END IF END IF END IF IF(LINEND.OR.ICB.EQ.0) THEN C End of the line or its part IF(IS.GT.KS) THEN IF( VV.GE.0.) THEN C Low-velocity EL1=EL1+AMAX1(0.,ELA-ELAOLD) EL2=EL2+AMAX1(0.,ELA-ELAOLD) AUX=AMIN1(PHI-PHIOLD,0.523599) EL2=EL2+0.5*ALOG(1.-SIN(2.*PHI0)*SIN(2.*AUX)) ELSE C High-velocity EL1=EL1 * +AMAX1(0.,ELA-ELAOLD+ALOG(ABS(COS(PHI-PHIOLD)))) EL2=EL2+AMAX1(0.,ELA-ELAOLD) * +ALOG(COS(AMIN1(PHI-PHIOLD,1.047198))) END IF IF(EL0.EQ.0..AND.VV.GT.0.) THEN ZON=ZON+1. END IF EL2=AMAX1(0.,EL2) * NUMSUM=NUMSUM+1 TTSUM =TTSUM +TT C TTSUM =TTSUM +S * ZONSUM=ZONSUM+ZON * ELASUM=ELASUM+ELA * PHISUM=PHISUM+PHI * EL1SUM=EL1SUM+EL1 EL2SUM=EL2SUM+EL2 EL0=EL0+AMAX1(0.,ELA-ELAOLD) IF(ABS(EL0-ELA).GT.0.000010) THEN C C MODLE2D-04 WRITE(TEXT,'(A,F8.3)') 'MODLE2D-04: Wrong EL0=',EL0 CALL ERROR(TEXT) END IF END IF END IF X1OLD=COOR1 X2OLD=COOR2 VOLD=VD(1) VVOLD=VV C *** End of numerical quadrature IF(LINEND) THEN C End of the line GO TO 71 END IF IF(ICB.EQ.0) THEN C Starting the new part of the line from the next point KS=IS+1 GO TO 20 END IF 70 CONTINUE C MODLE2D-05 CALL ERROR('MODLE2D-05: Line endpoint not reached') 71 CONTINUE 80 CONTINUE 85 CONTINUE C *** Results of numerical quadrature * AUX=TTSUM/FLOAT(NUMSUM) * ZON=ZONSUM/TTSUM * ELA=ELASUM/TTSUM * PHI=PHISUM/TTSUM * EL1=EL1SUM/TTSUM EL2=EL2SUM/TTSUM IF(FILE1.NE.' ') THEN IF(IA.EQ.0) THEN WRITE(LU1,'(9(A,F8.3))') '''LE'' /' * WRITE(LU1,'(2F9.3,A)') A,EL2,' M' ELSE * WRITE(LU1,'(2F9.3,A)') A,EL2,' T' END IF WRITE(LU1,'(9(A,F8.3))') ' ',A,' ',EL2,' /' END IF EL2MAX=AMAX1(EL2MAX,EL2) IF (FILE2.NE.' ') THEN READ(LU2,*) AW,W IF(ABS(AW-A).GT.ABS(0.01*DA)) THEN C C MODLE2D-06 CALL ERROR('MODLE2D-06: Wrong file with weights') C The angle AW read from the file MODLEW does not equal C to the angle A calculated from input parameters OA, DA. C Check the file MODLEW. END IF EL2=EL2*W ENDIF EL2MOD=EL2MOD+EL2 C *** WRITE(*,'(3(A,I3))') '+MODLE2D:',IA+1,' /',NA,' directions' 90 CONTINUE C C Closing input file MODLEW: IF (FILE2.NE.' ') CLOSE(LU2) C Closing output file with directional Lyapunov exponents: IF(FILE1.NE.' ') THEN WRITE(LU1,'(A)') '/' WRITE(LU1,'(A)') '/' CLOSE(LU1) END IF C C Output file with mean Lyapunov exponent: EL2MOD=EL2MOD/FLOAT(NA) CALL RSEP3T('MODLEM',FILE1,'modlem.out') IF(FILE1.NE.' ') THEN OPEN(LU1,FILE=FILE1) WRITE(LU1,'(A)') '/' AUX=EL2MOD+0.01 WRITE(LU1,'(9(A,F7.3))') '''',EL2MOD,''' ',OA, ' ',AUX ,' /' WRITE(LU1,'(9(A,F7.3))') ' ',OA ,' ',EL2MOD,' /' WRITE(LU1,'(9(A,F7.3))') ' ',OA+DA*FLOAT(NA-1),' ',EL2MOD,' /' WRITE(LU1,'(A)') '/' WRITE(LU1,'(A)') '/' CLOSE(LU1) END IF C C Output file with the frame: CALL RSEP3T('MODLEF',FILE1,'modlef.out') IF(FILE1.NE.' ') THEN CALL RSEP3R('LEMAX' ,EMAX ,EL2MAX) CALL RSEP3R('ALEMIN',ALEMIN,OA-0.5*DA) CALL RSEP3R('ALEMAX',ALEMAX,OA-0.5*DA+DA*FLOAT(NA)) OPEN(LU1,FILE=FILE1) WRITE(LU1,'(A)') '/' WRITE(LU1,'(9(A,F7.3))')'''',EL2MAX,''' ',OA,' ',EL2MAX-.03,' /' WRITE(LU1,'(9(A,F7.3))') ' ',ALEMIN,' ',EMAX,' /' WRITE(LU1,'(9(A,F7.3))') ' ',ALEMAX,' ',EMAX,' /' WRITE(LU1,'(9(A,F7.3))') ' ',ALEMAX,' ',0. ,' /' WRITE(LU1,'(9(A,F7.3))') ' ',ALEMIN,' ',0. ,' /' WRITE(LU1,'(9(A,F7.3))') ' ',ALEMIN,' ',EMAX,' /' WRITE(LU1,'(A)') '/' DO 92 I=0,INT(10.*EMAX+0.01) IF(MOD(I,10).EQ.0) THEN AUX=FLOAT(I/10) IF(I/10.LE.9) THEN WRITE(LU1,'(A,I1,9(A,F7.3))') * '''',I/10,''' ',ALEMIN-0.08,' ',AUX-0.03,' /' ELSE WRITE(LU1,'(A,I2,9(A,F7.3))') * '''',I/10,''' ',ALEMIN-0.16,' ',AUX-0.03,' /' END IF WRITE(LU1,'(9(A,F7.3))') ' ',ALEMIN ,' ',AUX ,' /' WRITE(LU1,'(9(A,F7.3))') ' ',ALEMIN+0.04,' ',AUX ,' /' WRITE(LU1,'(A)') '/' ELSE AUX=FLOAT(I)/10. WRITE(LU1,'(9(A,F7.3))') ''' '' ',ALEMIN-0.10,' ',AUX,' /' WRITE(LU1,'(9(A,F7.3))') ' ',ALEMIN ,' ',AUX,' /' WRITE(LU1,'(9(A,F7.3))') ' ',ALEMIN+0.02,' ',AUX,' /' WRITE(LU1,'(A)') '/' END IF 92 CONTINUE WRITE(LU1,'(A)') '/' CLOSE(LU1) END IF C WRITE(*,'(A)') '+MODLE2D: Done. ' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'length.for' C length.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