C
C *** OUTPUT OF LOVE WAVE PHASE AND GROUP VELOCITY AND 
C     AMPLITUDE SPECTRUM ON PLOTTER DATA FILE
C
      SUBROUTINE LVPTP(DOM,NOM,MOM,CFOM,NCFOM,CC,NSMD,MMD,UU,ALK,NADD)
C
C ** INPUT DATA **
C DOM     : DELTA OMEGA
C NOM     : TOTAL NUMBER OF CIRCULAR FREQUENCY
C MOM     : NOM IN MAIN PROGRAM DIMENSION
C         : IF MOM > 2**13(=8192), CAN NOT USE IN SUBROUTINE !!!
C ML      : NL IN MAIN 
C NL      : TOTAL NUMBER OF LAYERS           
C VS (I)  : VS OF Ith LAYER      
C MLE     : OUTPUT FINAL LAYER FOR ALK
C CFOM(I) : CUT-OFF OMEGA OF Ith MODE LOVE WAVE
C NCFOM   : NUMBER OF CUT-OFF OMEGA
C CC(I,J) : LOVE WAVE PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NSMD    : OUTPUT START LOVE MODE OF LOVE WAVE
C MMD     : NSMD IN MAIN PROGRAM DIMENSION
C UU(I,J) : GROUP VELOCITY OF Jth MODE LOVE WAVE AT Ith OMEGA
C ALK(I,J): AMPLITUDE SPECTRUM (AL/K) OF Jth MODE LOVE WAVE AT Ith OMEGA
C NADD    : =1: ADDING FROM FIRST TO FINAL MODE FOR ALK 
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD)
     *         ,UU(MOM,0:MMD),ALK(MOM,0:MMD)
      COMMON / MINMAX / VSMIN, VSMAX
C
      PAI2=8.D0*DATAN(1.D0)
C
C *** FINAL MODE FOR AMPLITUDE ***
C
      NAMPM=NCFOM
      IF(NADD.EQ.1)NAMPM=NSMD
C
C *** FIRST AND FINAL PERIOD, AND MAXIMUM AMPLITUDE SPECTRA AND PERIOD *
C
      FPR=PAI2/(DOM*DFLOAT(NOM))
      IF(NSMD.EQ.0)THEN
        EPR=PAI2/DOM
      ELSE
        EPR=PAI2/CFOM(NSMD)
      ENDIF
C
      ALKMAX=0.D0
      DO 210 MODE=NSMD,NAMPM
        DO 200 IOM=1,NOM
          ALKI=DABS(ALK(IOM,MODE))
          IF(ALKI.GT.ALKMAX)THEN
            ALKMAX=ALKI
            IMAX=IOM
          END IF
  200   CONTINUE
  210 CONTINUE
      PRMAX=PAI2/(DOM*DFLOAT(IMAX))
C
C *** OUTPUT OF DATA INFORMATION ***
C
      WRITE(*,*)
      WRITE(*,1000)
   10 WRITE(*,1010)
      WRITE(*,1020)DOM,NOM
      WRITE(*,1030)FPR,EPR
      WRITE(*,1040)VSMIN,VSMAX
      IF(ALKMAX.LT.0.001D0)THEN
      WRITE(*,1050)ALKMAX,PRMAX
      ELSE
      WRITE(*,1060)ALKMAX,PRMAX
      END IF
 1000 FORMAT(' From SUBROUTINE LVPTA : ')
 1010 FORMAT(/' DATA INFORMATION :')
 1020 FORMAT(' Delta Omega  = ',F8.4,' ,Omega Number = ',I5)
 1030 FORMAT(' First Period = ',F8.4,' ,Final Period = ',F8.4)
 1040 FORMAT(' Minimum Vs   = ',F8.3,' ,Maximum Vs   = ',F8.3)
 1050 FORMAT(' Maximum DABS(DAL) = ',D14.3, ' Period (sec) = ',F8.4)
 1060 FORMAT(' Maximum DABS(DAL) = ',F8.3, ' Period (sec) = ',F8.4)
C
C *** AXIS DATA INPUT ***
C
      WRITE(*,1200)
 1200 FORMAT(/' Axis Off ? ( Yes = 0 Enter ) :',\)
      READ(*,*)NAXY
C
      IF(NAXY.EQ.0)THEN
C
        OPEN(5,FILE='AXP.DAT')
        READ(5,1210)XOR,YOR,XFAC,YFAC,AMPFAC,XVAL,XLEN,XMIN
        CLOSE(5)
 1210   FORMAT(6D15.7)
C
      ELSE
C
        WRITE(*,1300)
        WRITE(*,1305)
        READ(*,*)XMIN
        WRITE(*,1310)
        READ(*,*)XVAL
        WRITE(*,1320)
        READ(*,*)XLEN
        XLEN=XLEN*100.0D0
        WRITE(*,1330)
        READ(*,*)NUX
        XFAC=XLEN/(XVAL-XMIN)
        WRITE(*,1335)
        READ(*,*)NXINT
        WRITE(*,1340)
        READ(*,*)YVAL
        WRITE(*,1350)
        READ(*,*)YLEN
        YLEN=YLEN*100.D0
        WRITE(*,1360)
        READ(*,*)NUY
        YFAC=YLEN/YVAL
        WRITE(*,1365)
        READ(*,*)NYVINT
        WRITE(*,1370)
        WRITE(*,1380)
        READ(*,*)YAVAL
        WRITE(*,1390)
        READ(*,*)YALEN
        YALEN=YALEN*100.D0
        WRITE(*,1400)
        READ(*,*)NAUY
        AMPFAC=YALEN/YAVAL
        WRITE(*,1405)
        READ(*,*)NYAINT
 1300   FORMAT(/' Axis Data Input for Phase and Group Velocity :')
 1305   FORMAT(/' Enter X-Axis Min (Min Period : sec)      :',\)
 1310   FORMAT( ' Enter X-Axis Max (Max Period : sec)      :',\)
 1320   FORMAT( ' Enter X-Axis Length (cm)                 :',\)
 1330   FORMAT( ' Enter X-Axis Division Number             :',\)
 1335   FORMAT( ' Enter Division Interval for X-Axis Scale :',\)
 1340   FORMAT(/' Enter Y-Axis Max Value (Max Vs : m/sec)  :',\)
 1350   FORMAT( ' Enter Y-Axis Length (cm)                 :',\)
 1360   FORMAT( ' Enter Y-Axis Division Number             :',\)
 1365   FORMAT( ' Enter Division Interval for Y-Axis Scale :',\)
 1370   FORMAT(/' Axis Data Input for Amplitude Spectrum :')
 1380   FORMAT(/' Enter Y-Axis Max Value (Amplitude factor):',\)
 1390   FORMAT( ' Enter Y-Axis Length (cm)                 :',\)
 1400   FORMAT( ' Enter Y-Axis Division Number             :',\)
 1405   FORMAT( ' Enter Division Interval for Y-Axis Scale :',\)
C
        WRITE(*,1500)
        READ(*,*)XOR
        XOR=(2.0D0+XOR)*100.0D0
        WRITE(*,1510)
        READ(*,*)YOR
        YOR=(2.0D0+YOR)*100.0D0
        WRITE(*,1520)
        READ(*,*)IYES
        IF(IYES.NE.0)GOTO 10
 1500   FORMAT(/' Enter X Origin Position ( cm ) :',\)
 1510   FORMAT( ' Enter Y Origin Position ( cm ) :',\)
 1520   FORMAT(/' All Input Data OK ? ( Yes = Enter 0 ) :',\)
C
      END IF
C
      OPEN(6,FILE='PS.DAT')
C
C *** OUTPUT OF PHASE AND GROUP VELOCITY PLOTTER DATA ***
C 
      IF(NAXY.EQ.0)GOTO 20
C
C   * PEN, TERM, AXIS AND CHARACTER *
C
      WRITE(6,2000)
      WRITE(6,2010)
C   * X-AXIS *
      WRITE(6,2020)XOR,YOR
      WRITE(6,2030)XLEN,NUX
      XINT=XLEN/DFLOAT(NUX)
      XINTV=(XVAL-XMIN)/DFLOAT(NUX)
      DO 90 INT=0,NUX,NXINT
        WRITE(6,2020)XOR+XINT*DFLOAT(INT)-100.0D0,YOR-60.0D0
        WRITE(6,2040)XMIN+XINTV*DFLOAT(INT)
   90 CONTINUE
      WRITE(6,2020)XOR+XLEN-270.0D0,YOR-120.0D0
      WRITE(6,2050)
C   * Y-AXIS (VEL) * 
      WRITE(6,2020)XOR,YOR
      WRITE(6,2060)YLEN,NUY
      YINT=YLEN/DFLOAT(NUY)
      YINTV=YVAL/DFLOAT(NUY)
      DO 95 INT=0,NUY,NYVINT
        WRITE(6,2020)XOR-40.0D0,YOR+YINT*DFLOAT(INT)-70.0D0
        WRITE(6,2070)YINTV*DFLOAT(INT)/1000.D0
   95 CONTINUE
      WRITE(6,2020)XOR-100.0D0,YOR+YLEN-200.0D0
      WRITE(6,2080)
C
 2000 FORMAT(2HJ1)
 2010 FORMAT(3H=<<)
 2020 FORMAT(1HM,F6.1,1H,F6.1,1H,)
 2030 FORMAT(3HX3,,F6.1,1H,I2,2H<<)
 2040 FORMAT(7HR0,S30,,1HP,F5.1,2H<<)
 2050 FORMAT(7HR0,S30,,14HPPERIOD(SEC)<<)
 2060 FORMAT(3HX2,,F6.1,1H,I2,2H<<)
 2070 FORMAT(9HR900,S30,,1HP,F4.1,2H<<)
 2080 FORMAT(9HR900,S30,,11HP(KM/SEC)<<)
C
   20 DO 100 MODE=NSMD,NCFOM
C
      CFPR=99999.D0
      IF(MODE.NE.0)CFPR=PAI2/CFOM(MODE)
C
C   * PHASE VELOCITY DATA PLOT *
C
      DO 105 IOM=NOM,1,-1
        PER=PAI2/(DOM*DFLOAT(IOM))
        IF(PER.GE.XMIN)GOTO 25
  105 CONTINUE
   25 FPER=PER
      IFOM=IOM
C
      WRITE(*,1600)MODE
      READ(*,*)NLIN
 1600 FORMAT(/' Line Type for Phase Velocity : Mode =',I5,
     *      /'   : Solid (0) or Dot (1) ? : ',\)
      WRITE(6,2094)
      IF(NLIN.EQ.1)WRITE(6,2095)
 2094 FORMAT(2HL0)
 2095 FORMAT(3HL3,,3HB30)
C
      WRITE(6,2100)XOR+(FPER-XMIN)*XFAC,YOR+CC(IFOM,MODE)*YFAC
      DO 110 IOM=IFOM,1,-1
        PER=PAI2/(DOM*DFLOAT(IOM))
        IF(PER.GE.CFPR)GO TO 30
        IF(PER.GE.XVAL)GO TO 30
        WRITE(6,2110)XOR+(PER-XMIN)*XFAC,YOR+CC(IOM,MODE)*YFAC
  110 CONTINUE
 2100 FORMAT(1HM,F6.1,1H,F6.1,1H,)
 2110 FORMAT(1HD,F6.1,1H,F6.1,1H,)
C
C   * GROUP VELOCITY DATA PLOT *
C
   30 WRITE(*,1700)MODE
      READ(*,*)NLIN
 1700 FORMAT(/' Line Type for Group Velocity : Mode =',I5,
     *      /'   : Solid (0) or Dot (1) ? : ',\)
      WRITE(6,2120)
      IF(NLIN.EQ.1)WRITE(6,2130)
 2120 FORMAT(2HL0)
 2130 FORMAT(3HL3,,3HB30)
C
      WRITE(6,2200)XOR+(FPER-XMIN)*XFAC,YOR+UU(IFOM,MODE)*YFAC
      DO 120 IOM=IFOM,1,-1
        PER=PAI2/(DOM*DFLOAT(IOM))
        IF(PER.GE.CFPR)GO TO 100
        IF(PER.GE.XVAL)GO TO 100
        WRITE(6,2110)XOR+(PER-XMIN)*XFAC,YOR+UU(IOM,MODE)*YFAC
  120 CONTINUE
 2200 FORMAT(1HM,F6.1,1H,F6.1,1H,)
 2210 FORMAT(1HD,F6.1,1H,F6.1,1H,)
C
  100 CONTINUE
C
C *** OUTPUT OF AMPLITUDE SPECTRUM PLOTTER DATA ***
C
      XORA=XOR+XLEN+300.D0
      IF(NAXY.EQ.0)GOTO 40
C
C   * PEN, TERM, AXIS AND CHARACTER *
C
      
      WRITE(6,2500)
      WRITE(6,2510)
C   * X-AXIS *
      WRITE(6,2520)XORA,YOR
      WRITE(6,2530)XLEN,NUX
      XINT=XLEN/DFLOAT(NUX)
      XINTV=(XVAL-XMIN)/DFLOAT(NUX)
      DO 97 INT=0,NUX,NXINT
        WRITE(6,2520)XORA+XINT*DFLOAT(INT)-100.0D0,YOR-60.0D0
        WRITE(6,2540)XMIN+XINTV*DFLOAT(INT)
   97 CONTINUE
      WRITE(6,2520)XORA+XLEN-270.0D0,YOR-120.0D0
      WRITE(6,2550)
C   * Y-AXIS (AMP) *
      WRITE(6,2520)XORA,YOR
      WRITE(6,2560)YALEN,NAUY
      YINT=YALEN/DFLOAT(NAUY)
      YINTV=YAVAL/DFLOAT(NAUY)
      DO 98 INT=0,NAUY,NYAINT
        WRITE(6,2520)XORA-40.0D0,YOR+YINT*DFLOAT(INT)-70.0D0
        YAMP=YINTV*DFLOAT(INT)
C        IF(YAMP
        WRITE(6,2570)YAMP
   98 CONTINUE
      WRITE(6,2520)XORA-100.0D0,YOR+YALEN-480.0D0
      WRITE(6,2580)
C
 2500 FORMAT(2HJ1)
 2510 FORMAT(3H=<<)
 2520 FORMAT(1HM,F6.1,1H,F6.1,1H,)
 2530 FORMAT(3HX3,,F6.1,1H,I2,2H<<)
 2540 FORMAT(7HR0,S30,,1HP,F5.1,2H<<)
 2550 FORMAT(7HR0,S30,,14HPPERIOD(SEC)<<)
 2560 FORMAT(3HX2,,F6.1,1H,I2,2H<<)
 2570 FORMAT(9HR900,S30,,1HP,F4.1,2H<<)
 2575 FORMAT(9HR900,S30,,1HP,D11.3,2H<<)
 2580 FORMAT(9HR900,S25,,4HQ25,,23HPAMPLIFICATION FACTOR<<,3HQ30)
C
   40 DO 300 MODE=NSMD,NAMPM
C
      CFPR=99999.D0
      IF(MODE.NE.0)CFPR=PAI2/CFOM(MODE)
C
C   * AMPLITUDE SPECTRUM DATA PLOT *
C
      WRITE(*,1800)MODE
      READ(*,*)NLIN
 1800 FORMAT(/' Line Type for Amplitude Spectrum : Mode =',I5,
     *      /'   : Solid (0) or Dot (1) ? : ',\)
      WRITE(6,2600)
      IF(NLIN.EQ.1)WRITE(6,2610)
 2600 FORMAT(2HL0)
 2610 FORMAT(3HL3,,3HB30)
C
      WRITE(6,2700)XORA+(FPER-XMIN)*XFAC
     *              ,YOR+DABS(ALK(IFOM,MODE)*AMPFAC)
C
      DO 310 IOM=IFOM,1,-1
        PER=PAI2/(DOM*DFLOAT(IOM))
        IF(PER.GE.CFPR)GO TO 300
        IF(PER.GE.XVAL)GO TO 300
        WRITE(6,2710)XORA+(PER-XMIN)*XFAC,YOR+DABS(ALK(IOM,MODE)*AMPFAC)
  310 CONTINUE
 2700 FORMAT(1HM,F6.1,1H,F6.1,1H,)
 2710 FORMAT(1HD,F6.1,1H,F6.1,1H,)
C
  300 CONTINUE
C
      WRITE(6,2800)XORA+PER*XFAC+500.D0,YOR
 2800 FORMAT(1HM,F6.1,1H,F6.1,1H,)
C
      CLOSE(6)
C
C *** AXIS DATA OUTPUT TO FILE ***
C
      IF(NAXY.NE.0)THEN
        OPEN(5,FILE='AXP.DAT')
        WRITE(5,3000)XOR,YOR,XFAC,YFAC,AMPFAC,XVAL,XLEN,XMIN
        CLOSE(5)
 3000   FORMAT(6D15.7)
      END IF
C
      WRITE(*,3200)
 3200 FORMAT(/' USE FOR PLOTTER : A> COPY PS.DAT PRN')
      WRITE(*,*)
C
      RETURN
      END
