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.50
C Date: 2001, June 1
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 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.
C             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 avaraging
C             directional Lyapunov exponents with a unit weight.
C             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
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
      CHARACTER*30 TEXT
      PARAMETER (LU1=1)
      LOGICAL LINEND
      REAL COOR(3),H(3)
      REAL UP(10),US(10),VD(10),AUX0,AUX1,AUX2,AUX3,AUX4
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     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
        EL2MOD=EL2MOD+EL2
        EL2MAX=AMAX1(EL2MAX,EL2)
C       ***
        WRITE(*,'(3(A,I3))') '+MODLE2D:',IA+1,' /',NA,' directions'
   90 CONTINUE
C
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