CC *** READ OF SOIL DATA FROM DATA FILE ***
C
      SUBROUTINE RINPUT(DOM,NOM,ML,NL,DNS,VP,VS,DEP,NDOM,NSMD,NEMD,NDC)
C
C ** INPUT DATA **
C DOM   : INCREMENT OF CIRCULAR FREQUENCY (DELTA OMEGA)
C NOM   : TOTAL NUMBER OF CIRCULAR FREQUENCIES
C ML    : NUMBER OF SOIL LAYERS IN MAIN DIMENSION
C NL    : NUMBER OF SOIL LAYERS
C DNS(I): DENSITY OF Ith LAYER 
C VP (I): VP OF Ith LAYER                              
C VS (I): VS OF Ith LAYER                              
C DEP(I): DEPTH (THICKNESS) OF Ith LAYER             
C NDOM  : DIVISION NUMBER BETWEEN 0 AND MAX OMEGAS FOR CUT-OFF OMEGAS
C NSMD  : FIRST RAYLEIGH MODE TO CALCURATE
C NEMD  : FINAL RAYLEIGH MODE TO CALCULATE
C NDC   : DIVISION NUMBER BETWEEN MIN AND MAX VS FOR PHASE VELOCITIES
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*30
      DIMENSION DNS(ML),VP(ML),VS(ML),DEP(ML)
      COMMON / MINMAX / VSMIN,VSMAX,CCMIN,CCMAX,MINIL,MAXIL
C
C *** SOIL DATA INPUT FROM FILE DATA ***
C
      WRITE(*,*)
      WRITE(*,1000)
   10 WRITE(*,1010)
      READ(*,'(A)',ERR=10)FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 10
 1000 FORMAT( ' From SUBROUTINE RINPUT :')
 1010 FORMAT(/' Enter Soil Data File Name ',\)

      OPEN (7,FILE=FNAME)
      READ(7,1200)NL
      DO 100 IL=1,NL
        READ(7,1210)DNS(IL),VP(IL),VS(IL),DEP(IL),QF
  100 CONTINUE
      CLOSE(7)
 1200 FORMAT(I10)
 1210 FORMAT(5F10.0)
C
C *** DATA ECHO ***
C
      WRITE(*,2000)NL
      DO 200 IL=1,NL
        WRITE(*,2010)IL,DNS(IL),VP(IL),VS(IL),DEP(IL)
  200 CONTINUE
 2000 FORMAT(/' Total Layers Number            =',I5)
 2010 FORMAT( ' Layer Num., Density, Vp, Vs, Depth =',I3,4F10.3)
C
C *** MIN AND MAX VS AND SOIL DATA CHECK ***
C
      VSMIN=VS(NL)
      VSMAX=0.0D0
      VPMAX=0.0D0 
      MINIL=1
      MAXIL=1
      DO 300 IL=1,NL
        IF(VS(IL).LE.VSMIN)THEN
          VSMIN=VS(IL)
          MINIL=IL
        END IF
        IF(VS(IL).GE.VSMAX)THEN
          VSMAX=VS(IL)
          MAXIL=IL
        END IF
        VPMAX=MAX(VPMAX,VP(IL))
  300 CONTINUE
C
      IF(VSMAX.NE.VS(NL))THEN
        WRITE(*,*)' FROM SUBROUTINE RINPUT :'
        WRITE(*,*)' VS(NL).NE.VSMAX, THUS STOP !!!'
        STOP
      END IF
      IF(VPMAX.NE.VP(NL))THEN
        WRITE(*,*)' FROM SUBROUTINE RINPUT :'
        WRITE(*,*)' VP(NL).NE.VPMAX, THUS STOP !!!'
        STOP
      END IF
C
C *** MINUMUM AND MAXIMUM PHASE VELOCITIES ***
C
      CCMIN=FNCPH1(VP(MINIL),VS(MINIL))
      CCMAX=FNCPH1(VP(MAXIL),VS(MAXIL))
C
C *** DATA INPUT FROM SCREEN ***
C
   20 WRITE(*,3000)
      READ(*,*)DOM,NOM
      OMMAX=DOM*DFLOAT(NOM)
      WRITE(*,3010)0.D0,OMMAX
C
      WRITE(*,3020)
      READ(*,*)NSMD
      WRITE(*,3030)
      READ(*,*)NEMD
      IF(NEMD.EQ.0)GOTO 30
C
      WRITE(*,3040)
      READ(*,*)NDOM
C
   30 WRITE(*,3050)CCMIN,VSMAX
      WRITE(*,3060)
      READ(*,*)NDC
C
      WRITE(*,3070)
      READ(*,*)IYES
      IF(IYES.NE.0)GOTO 20
C
 3000 FORMAT(/' Enter Delta Omega and Omega Number         ? :',\)
 3010 FORMAT( ' Minimum Omega and Maximum Omega = ',2F8.3)
 3020 FORMAT(/' Enter First Rayleigh Mode (>= 0  )         ? :',\)
 3030 FORMAT( ' Enter Final Rayleigh Mode                  ? :',\)
 3040 FORMAT(/' To Get Cut-Off Omegas for higher mode,  Enter',
     *       /' Division Number between Min and Max Omegas ? :',\)
 3050 FORMAT(/' Minimum C and Maximum Vs = ',2F10.3)
 3060 FORMAT( ' To Get Phase Velocities (C), Enter',
     *       /' Division Number between Min and Max Cs     ? :',\)
 3070 FORMAT(/' Input Data OK ? (yes = Enter 0)            ? :',\)
      WRITE(*,*)
C
      RETURN
      END
C
C
C
C *** CAL. OF RAYLEIGH PHASE WAVE VELOCITY IN A MEDIUM OF SINGLE LAYER ***
C
      REAL FUNCTION FNCPH1*8(VP,VS)
      IMPLICIT REAL*8(A-H,O-Z)
C
      DZ=0.0001D0
C
      SP=VS/VP
      IF(SP.GE.0.0D0) Z= 0.91263D0
      IF(SP.GE.0.2D0) Z= 0.90767D0
      IF(SP.GE.0.4D0) Z= 0.88899D0
      IF(SP.GE.0.6D0) Z= 0.83575D0
      IF(SP.GE.0.7D0) Z= 0.77062D0
      IF(SP.GE.0.8D0) Z= 0.63338D0
C
      F0=FNCFR1(SP,Z)
  100 Z=Z-DZ
      F1=FNCFR1(SP,Z)
      AF=F0/F1
      F00=F0
      F0 =F1
      IF(AF.GE.0.D0)GOTO 100
      FZ=Z+DZ*F0/(F0-F00)
C
      FNCPH1=VS*DSQRT(FZ)
C
      RETURN
      END
C
C
      REAL FUNCTION FNCFR1*8(SP,Z)
      IMPLICIT REAL*8(A-H,O-Z)
C
      Z2=Z*Z
      Z3=Z2*Z
C
      F1=16.D0*(Z-1.D0)*SP*SP
      FNCFR1=F1-Z3+8.D0*Z2-24.D0*Z+16.D0
C
      RETURN
      END
C
C
C
C *** OUTPUT OF PHASE & GROUP VELOCITIES AND 
C     AMPLITUDE SPECTRA ON SCREEN
C
      SUBROUTINE ROUT(DOM,NOM,MOM,ML,NL,VS
     *               ,CFOM,NCFOM,CC,NSMD,MMD,UU,ARV,ARH,UW0)
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 CFOM(I) : CUT-OFF OMEGA OF Ith MODE RAYLEIGH WAVE
C NCFOM   : NUMBER OF CUT-OFF OMEGA
C CC(I,J) : RAYLEIGH WAVE PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NSMD    : OUTPUT START RAYLEIGH MODE
C MMD     : NMD IN MAIN PROGRAM DIMENSION
C UU(I,J) : GROUP VELOCITY OF Jth MODE RAYLEIGH WAVE AT Ith OMEGA
C ARV(I,J): AMPLITUDE SPECTRUM (DISPLACEMENT) DUE TO VERTICAL UNIT POINT-FOURCE
c           OF Jth MODE AT Ith OMEGA (REAL)
C ARH(I,J): AMPLITUDE SPECTRUM (DISPLACEMENT) DUE TO HORIZONTAL 
C           UNIT POINT-FOURCE OF Jth MODE AT Ith OMEGA (REAL)
C UW0(I,J): AMPLITUDE RATIO ( IMAG(U0/W0 ) OF Jth MODE RAYLEIGH WAVE 
C           AT Ith OMEGA
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD)
     *         ,UU(MOM,0:MMD),ARV(MOM,0:MMD),ARH(MOM,0:MMD)
     *         ,UW0(MOM,0:MMD),VS(ML)
      COMMON / MINMAX / VSMIN,VSMAX,CCMIN,CCMAX,MINIL,MAXIL
C
      WRITE(*,1000)
      WRITE(*,1010)NCFOM
 1000 FORMAT(/' From SUBROUTINE ROUT :')
 1010 FORMAT(/' Highest Mode =',I5)
C
      DO 100 MODE=NSMD,NCFOM
C
        CFCC=VS(NL)
        IF(MODE.EQ.0)CFCC=CCMAX
        WRITE(*,1100)MODE
        WRITE(*,1110)CFOM(MODE),CFCC
        WRITE(*,1120)
C
        ICN=1
        DO 120 IOM=1,NOM
C
         OM=DOM*FLOAT(IOM)
         IF(OM.LT.CFOM(MODE))GOTO 120
         PER=6.2831853070D0/OM
         WRITE(*,1130)ICN,OM,PER,CC(IOM,MODE),UU(IOM,MODE)
     *          ,ARV(IOM,MODE),ARH(IOM,MODE),UW0(IOM,MODE)
         ICN=ICN+1
C
  120   CONTINUE
C
  100 CONTINUE
C
 1100 FORMAT(//' Rayleigh Mode =',I5)
 1110 FORMAT(' Cut-Off Omega =',F10.5,',  Phase and Group Vel. =',F9.3/)
 1120 FORMAT(' Num,Omega,Period, Phase V., Group V.,  Ver.Spec'
     *       ,',  Hor.Spec,  Imag(U/W)'/)
 1130 FORMAT(I3,F7.3,F7.2,F10.2,F10.2,1X,D10.3,1X,D10.3,D11.3)  
C
      WRITE(*,*)
C
      RETURN
      END
C
C
C
C *** PLOT OF PHASE AND GROUP VELOCITIES AND 
C     AMPLITUDE SPECTUM DUE TO VERTICAL FOURCE ***
C
      SUBROUTINE RPLOTV(DOM,NOM,MOM,ML,NL,VS
     *                 ,CFOM,NCFOM,CC,NSMD,MMD,UU,ARV,UW0)
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 CFOM(I) : CUT-OFF OMEGA OF Ith MODE RAYLEIGH WAVE
C NCFOM   : NUMBER OF CUT-OFF OMEGA
C CC(I,J) : RAYLEIGH WAVE PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NSMD    : OUTPUT START RAYLEIGH MODE
C MMD     : NSMD IN MAIN PROGRAM DIMENSION
C UU(I,J) : GROUP VELOCITY OF Jth MODE RAYLEIGH WAVE AT Ith OMEGA
C ARV(I,J): AMPLITUDE SPECTRUM (DISPLACEMENT) OF Jth MODE AT Ith OMEGA
C           AT Ith OMEGA
C
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*2 NP
      DIMENSION CFOM(0:MMD),CC(MOM,0:MMD)
     *         ,UU(MOM,0:MMD),ARV(MOM,0:MMD),UW0(MOM,0:MMD),VS(ML)
      COMMON / MINMAX / VSMIN,VSMAX,CCMIN,CCMAX,MINIL,MAXIL
C
C *** FINAL OMEGA ***
C
      EOM=DOM*DFLOAT(NOM)
C
C *** OUTPUT OF DATA INFORMATION ***
C
      WRITE(*,*)
      WRITE(*,1000)
   10 WRITE(*,1010)
      WRITE(*,1020)DOM,NOM
      WRITE(*,1030)CFOM(NSMD),EOM
      WRITE(*,1040)CCMIN,VSMAX
 1000 FORMAT(' From SUBROUTINE RPLOTV : ')
 1010 FORMAT(/' DATA INFORMATION :')
 1020 FORMAT(' Delta Omega = ',F8.4,' ,Omega Number = ',I5)
 1030 FORMAT(' First Omega = ',F8.4,' ,Final Omega  = ',F8.4)
 1040 FORMAT(' Minimum C   = ',F8.3,' ,Maximum Vs   = ',F8.3)
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='AXV.DAT')
          READ(5,1210)XOR,YOR,XFAC,YFAC,AMPFAC,XVAL,IYAMP
        CLOSE(5)
 1210   FORMAT(6D15.7,I5)
C
      ELSE
C
        WRITE(*,1300)
        WRITE(*,1310)
        READ(*,*)XVAL
        WRITE(*,1320)
        READ(*,*)XLEN
        XLEN=XLEN*100.0D0
        WRITE(*,1330)
        READ(*,*)NUX
        XFAC=XLEN/XVAL
        WRITE(*,1340)
        READ(*,*)YVAL
        WRITE(*,1350)
        READ(*,*)YLEN
        YLEN=YLEN*100.D0
        WRITE(*,1360)
        READ(*,*)NUY
        YFAC=YLEN/YVAL
        WRITE(*,1370)
        WRITE(*,1380)
        READ(*,*)YAMP
        YAMP=YAMP*100.D0
        WRITE(*,1390)
        READ(*,*)IYAMP
 1300   FORMAT(/' Axis Data Input for Phase and Group Velocity :')
 1310   FORMAT(/' Enter X-Axis Max (Max Omega: rad/sec)   :',\)
 1320   FORMAT( ' Enter X-Axis Length (cm)                :',\)
 1330   FORMAT( ' Enter X-Axis Division Number            :',\)
 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            :',\)
 1370   FORMAT(/' Axis Data Input for Amplitude Spectrum :')
 1380   FORMAT(/' Enter Y-Axis Length (cm)                :',\)
 1390   FORMAT(/' Superpose Amp W and U ? (Yes = 0 Enter) :',\)
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
C *** PLOTTER DATA FILE OPEN ***
C
      OPEN(6,FILE='PV.DAT')
C
C    ** OUTPUT OF PHASE AND GROUP VELOCITY PLOTTER DATA ***
C 
        WRITE(6,2000)
C
        IF(NAXY.EQ.0)GOTO 20
C
C     * TERM, AXIS AND CHARACTER *
C
        WRITE(6,2010)
        WRITE(6,2020)XOR,YOR
        WRITE(6,2030)XLEN,NUX
        WRITE(6,2020)XOR-110.0D0,YOR-60.0D0
        WRITE(6,2040)0.0D0
        WRITE(6,2020)XOR+XLEN-100.0D0,YOR-60.0D0
        WRITE(6,2040)XVAL
        WRITE(6,2020)XOR+XLEN-220.0D0,YOR-120.0D0
        WRITE(6,2050)
        WRITE(6,2020)XOR,YOR
        WRITE(6,2060)YLEN,NUY
        WRITE(6,2020)XOR-40.0D0,YOR-80.0D0
        WRITE(6,2070)0.0D0
        WRITE(6,2020)XOR-40.0D0,YOR+YLEN-80.0D0
        WRITE(6,2080)YVAL/1000.D0
        WRITE(6,2020)XOR-100.0D0,YOR+YLEN-200.0D0
        WRITE(6,2090)
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,,12HP(RAD/SEC)<<)
 2060   FORMAT(3HX2,,F6.1,1H,I2,2H<<)
 2070   FORMAT(9HR900,S30,,1HP,F4.1,2H<<)
 2080   FORMAT(9HR900,S30,,1HP,F4.1,2H<<)
 2090   FORMAT(9HR900,S30,,11HP(KM/SEC)<<)
C
   20   DO 100 MODE=NSMD,NCFOM
C
C     * PHASE VELOCITY DATA PLOT *
C
          WRITE(*,2001)MODE
          READ(*,*)NLIN
 2001     FORMAT(/' Line Type for Phase Velocity : Mode =',I5,
     *           /'   : Solid (0) or Dot (1) ? :',\)
          WRITE(6,2002)
          IF(NLIN.EQ.1)WRITE(6,2003)
 2002     FORMAT(2HL0)
 2003     FORMAT(3HL3,,3HB30)
          WRITE(*,2004)
          READ(*,'(A)')NP
 2004     FORMAT(/' Pen Number (J1-J9) ? :',\)
          WRITE(6,2005)NP
 2005     FORMAT(A)
C
          YV=VS(NL)
          IF(MODE.EQ.0)YV=CCMAX
C
          WRITE(6,2100)XOR+CFOM(MODE)*XFAC,YOR+YV*YFAC
          DO 110 IOM=1,NOM
            OMX=DOM*DFLOAT(IOM)
            IF(OMX.LT.CFOM(MODE))GOTO 110
            WRITE(6,2110)XOR+OMX*XFAC,YOR+CC(IOM,MODE)*YFAC
            IF(OMX.GE.XVAL)GO TO 30
  110     CONTINUE
C
 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(*,2101)MODE
          READ(*,*)NLIN
 2101     FORMAT(/' Line Type for Group Velocity : Mode =',I5,
     *           /'   : Solid (0) or Dot (1) ? :',\)
          WRITE(6,2102)
          IF(NLIN.EQ.1)WRITE(6,2103)
 2102     FORMAT(2HL0)
 2103     FORMAT(3HL3,,3HB30)
          WRITE(*,2104)
          READ(*,'(A)')NP
 2104     FORMAT(/' Pen Number (J1-J9) ? :',\)
          WRITE(6,2105)NP
 2105     FORMAT(A)
C
          WRITE(6,2200)XOR+CFOM(MODE)*XFAC,YOR+YV*YFAC
          DO 120 IOM=1,NOM
            OMX=DOM*DFLOAT(IOM)
            IF(OMX.LT.CFOM(MODE))GOTO 120
            WRITE(6,2110)XOR+OMX*XFAC,YOR+UU(IOM,MODE)*YFAC
            IF(OMX.GE.XVAL)GO TO 100
  120     CONTINUE
C
 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 W SPECTRA PLOTTER DATA **
C
        XORA=XOR+XVAL*XFAC+200.D0
        IF(NAXY.EQ.0)GOTO 40
C
C     * MAXIMUM AMPLITUDE SPECTRA *
C
        IOMMAX=0
        MDMAX=0
        ARVMAX=0.D0
        DO 210 MODE=NSMD,NCFOM
          DO 200 IOM=1,NOM
            ARVI=DABS(ARV(IOM,MODE))
            IF(ARVI.GT.ARVMAX)THEN
              IOMMAX=IOM
              MDMAX=MODE
              ARVMAX=ARVI
            END IF
            ARVI=DABS(ARV(IOM,MODE)*UW0(IOM,MODE))
            IF(ARVI.GT.ARVMAX)THEN
              IOMMAX=IOM
              MDMAX=MODE
              ARVMAX=ARVI
            END IF
  200     CONTINUE
  210   CONTINUE
        AMPFAC=YAMP/ARVMAX
C
C     * PEN, TERM, AXIS AND CHARACTER *
C
        WRITE(6,2500)
        WRITE(6,2510)
        WRITE(6,2520)XORA,YOR
        WRITE(6,2530)XLEN,NUX
        WRITE(6,2520)XORA-110.0D0,YOR-60.0D0
        WRITE(6,2540)0.0D0
        WRITE(6,2520)XORA+XLEN-100.0D0,YOR-60.0D0
        WRITE(6,2540)XVAL
        WRITE(6,2520)XORA+XLEN-220.0D0,YOR-120.0D0
        WRITE(6,2550)
        WRITE(6,2520)XORA,YOR
        WRITE(6,2560)YAMP,1
        WRITE(6,2520)XORA-40.0D0,YOR-80.0D0
        WRITE(6,2570)0.0D0
        WRITE(6,2520)XORA-40.0D0,YOR+YAMP-250.0D0
        WRITE(6,2580)ARVMAX
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,,12HP(RAD/SEC)<<)
 2560   FORMAT(3HX2,,F6.1,1H,I2,2H<<)
 2570   FORMAT(9HR900,S30,,1HP,F4.1,2H<<)
 2580   FORMAT(9HR900,S30,,1HP,D11.3,2H<<)
C
C     * AMPLITUDE SPECTRA DATA PLOT ( W ) *
C
   40   DO 300 MODE=NSMD,NCFOM
C
          WRITE(*,2501)MODE
          READ(*,*)NLIN
 2501     FORMAT(/' Line Type for Vertical Amp.Spec.(ARV) : Mode =',I5,
     *           /'   : Solid (0) or Dot (1) ? :',\)
          WRITE(6,2502)
          IF(NLIN.EQ.1)WRITE(6,2503)
 2502     FORMAT(2HL0)
 2503     FORMAT(3HL3,,3HB30)
          WRITE(*,2504)
          READ(*,'(A)')NP
 2504     FORMAT(/' Pen Number (J1-J9) ? :',\)
          WRITE(6,2505)NP
 2505     FORMAT(A)
C
          WRITE(6,2700)XORA+CFOM(MODE)*XFAC,YOR
          DO 310 IOM=1,NOM
            OMX=DOM*DFLOAT(IOM)
            IF(OMX.LT.CFOM(MODE))GOTO 310
            WRITE(6,2710)XORA+OMX*XFAC,YOR+DABS(ARV(IOM,MODE)*AMPFAC)
            IF(OMX.GE.XVAL)GO TO 300
  310     CONTINUE
C
 2700     FORMAT(1HM,F6.1,1H,F6.1,1H,)
 2710     FORMAT(1HD,F6.1,1H,F6.1,1H,)
C
  300   CONTINUE
C
C    ** OUTPUT OF AMPLITUDE U SPECTRA PLOTTER DATA **
C
        IF(IYAMP.EQ.0)THEN
          WRITE(6,2800)
          GOTO 50
        END IF
 2800   FORMAT(2HJ2)
C
        XORA=XORA+XVAL*XFAC+200.D0
        IF(NAXY.EQ.0)GOTO 50
C
C     * PEN, TERM, AXIS AND CHARACTER *
C
      
        WRITE(6,3000)
        WRITE(6,3010)
        WRITE(6,3020)XORA,YOR
        WRITE(6,3030)XLEN,NUX
        WRITE(6,3020)XORA-110.0D0,YOR-60.0D0
        WRITE(6,3040)0.0D0
        WRITE(6,3020)XORA+XLEN-100.0D0,YOR-60.0D0
        WRITE(6,3040)XVAL
        WRITE(6,3020)XORA+XLEN-220.0D0,YOR-120.0D0
        WRITE(6,3050)
        WRITE(6,3020)XORA,YOR
        WRITE(6,3060)YAMP,1
        WRITE(6,3020)XORA-40.0D0,YOR-80.0D0
        WRITE(6,3070)0.0D0
        WRITE(6,3020)XORA-40.0D0,YOR+YAMP-250.0D0
        WRITE(6,3080)ARVMAX
C
 3000   FORMAT(2HJ1)
 3010   FORMAT(3H=<<)
 3020   FORMAT(1HM,F6.1,1H,F6.1,1H,)
 3030   FORMAT(3HX3,,F6.1,1H,I2,2H<<)
 3040   FORMAT(7HR0,S30,,1HP,F5.1,2H<<)
 3050   FORMAT(7HR0,S30,,12HP(RAD/SEC)<<)
 3060   FORMAT(3HX2,,F6.1,1H,I2,2H<<)
 3070   FORMAT(9HR900,S30,,1HP,F4.1,2H<<)
 3080   FORMAT(9HR900,S30,,1HP,D11.3,2H<<)
C
C     * AMPLITUDE SPECTRA DATA PLOT ( U ) *
C
   50   DO 500 MODE=NSMD,NCFOM
C
          WRITE(*,3001)MODE
          READ(*,*)NLIN
 3001     FORMAT(/' Line Type for Horizontal Amp.Spec.(ARH) : Mode =',I5,
     *           /'   : Solid (0) or Dot (1) ? :',\)
          WRITE(6,3002)
          IF(NLIN.EQ.1)WRITE(6,3003)
 3002     FORMAT(2HL0)
 3003     FORMAT(3HL3,,3HB30)
          WRITE(*,3004)
          READ(*,'(A)')NP
 3004     FORMAT(/' Pen Number (J1-J9) ? :',\)
          WRITE(6,3005)NP
 3005     FORMAT(A)
C
          WRITE(6,3210)XORA+CFOM(MODE)*XFAC,YOR
          DO 510 IOM=1,NOM
            OMX=DOM*DFLOAT(IOM)
            IF(OMX.LT.CFOM(MODE))GOTO 510
            ARVU=DABS(ARV(IOM,MODE)*AMPFAC*UW0(IOM,MODE))
            WRITE(6,3220)XORA+OMX*XFAC,YOR+ARVU
            IF(OMX.GE.XVAL)GO TO 500
  510     CONTINUE
C
 3210     FORMAT(1HM,F6.1,1H,F6.1,1H,)
 3220     FORMAT(1HD,F6.1,1H,F6.1,1H,)
C
  500   CONTINUE
C
C     * SKIP OF PLOTTER PEN *
C
        WRITE(6,3300)XORA+XVAL*XFAC+500.D0,YOR
 3300   FORMAT(1HM,F6.1,1H,F6.1,1H,)
C
C *** POTTER DATA FILE CLOSE ***
C
      CLOSE(6)
C
C *** AXIS DATA OUTPUT TO FILE ***
C
      IF(NAXY.NE.0)THEN
        OPEN(5,FILE='AXV.DAT')
        WRITE(5,3500)XOR,YOR,XFAC,YFAC,AMPFAC
     *              ,XVAL,IYAMP,XLEN,YLEN,YAMP,NUX
        CLOSE(5)
 3500   FORMAT(6D15.7,I5,3D15.7,I5)
      END IF
C
      WRITE(*,3700)
 3700 FORMAT(/' USE FOR PLOTTER : A> COPY PV.DAT PRN')
      WRITE(*,*)
C
      RETURN
      END
C
C
C
C *** PLOT OF AMPLITUDE SPECTRUM DUE TO HORIZONTAL FORCE ***
C
      SUBROUTINE RPLOTH(DOM,NOM,MOM,CFOM,NCFOM,NSMD,MMD,ARH,UW0)
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 CFOM(I) : CUT-OFF OMEGA OF Ith MODE RAYLEIGH WAVE
C NCFOM   : NUMBER OF CUT-OFF OMEGA
C NSMD    : OUTPUT START RAYLEIGH MODE
C MMD     : NSMD IN MAIN PROGRAM DIMENSION
C ARH(I,J): AMPLITUDE SPECTRUM (DISPLACEMENT) OF Jth MODE AT Ith OMEGA
C UW0(I,J): AMPLITUDE RATIO (IMAG(RK/RL)) OF Jth MODE AT Ith OMEGA
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION CFOM(0:MMD),ARH(MOM,0:MMD),UW0(MOM,0:MMD)
      COMMON / MINMAX / VSMIN,VSMAX,CCMIN,CCMAX,MINIL,MAXIL
C
C *** AXIS DATA INPUT ***
C
      WRITE(*,1200)
      WRITE(*,1201)
 1200 FORMAT(/' From SUBROUTINE RPLOTH :')
 1201 FORMAT(/' Axis Off ? (Second Time ? ; Yes = 0) :',\)
      READ(*,*)NAXY
C
      IF(NAXY.EQ.0)THEN
C
        OPEN(5,FILE='AXH.DAT')
          READ(5,1210)XOR,YOR,XFAC,YFAC,AMPFAC,XVAL,IYAMP,YLEN,NVH
        CLOSE(5)
 1210   FORMAT(6D15.7,I5,D15.7,I5)
C
      ELSE
C
        WRITE(*,1202)
 1202   FORMAT(/' Superpose ARV and ARH ? (Yes = 0) :',\)
        READ(*,*)NVH
C
        OPEN(5,FILE='AXV.DAT')
          READ(5,1211)XOR,YOR,XFAC,YFAC,AMPFAC
     *               ,XVAL,IYAMP,XLEN,YLEN,YAMP,NUX
        CLOSE(5)
 1211   FORMAT(6D15.7,I5,3D15.7,I5)
C
      END IF
C
C *** PLOTTER DATA FILE OPEN ***
C
      OPEN(6,FILE='PH.DAT')
C
C    ** OUTPUT OF AMPLITUDE W SPECTRA PLOTTER DATA **
C
        XORA=XOR+XVAL*XFAC+200.D0
        YORO=YOR
        IF(NVH.NE.0)YOR=YOR+YLEN+200.D0
        IF((NAXY.NE.0).AND.(NVH.NE.0))THEN
C
C       * MAXIMUM AMPLITUDE SPECTRA *
C
          IOMMAX=0
          MDMAX=0
          ARHMAX=0.D0
          DO 210 MODE=NSMD,NCFOM
            DO 200 IOM=1,NOM
              ARHI=DABS(ARH(IOM,MODE))
              IF(ARHI.GT.ARHMAX)THEN
                IOMMAX=IOM
                MDMAX=MODE
                ARHMAX=ARHI
              END IF
              ARHI=DABS(ARH(IOM,MODE)*UW0(IOM,MODE))
              IF(ARHI.GT.ARHMAX)THEN
                IOMMAX=IOM
                MDMAX=MODE
                ARHMAX=ARHI
              END IF
  200       CONTINUE
  210     CONTINUE
          AMPFAC=YAMP/ARHMAX
C
C       * PEN, TERM, AXIS AND CHARACTER *
C
          WRITE(6,2500)
          WRITE(6,2510)
          WRITE(6,2520)XORA,YOR
          WRITE(6,2530)XLEN,NUX
          WRITE(6,2520)XORA-110.0D0,YOR-60.0D0
          WRITE(6,2540)0.0D0
          WRITE(6,2520)XORA+XLEN-100.0D0,YOR-60.0D0
          WRITE(6,2540)XVAL
          WRITE(6,2520)XORA+XLEN-220.0D0,YOR-120.0D0
          WRITE(6,2550)
          WRITE(6,2520)XORA,YOR
          WRITE(6,2560)YAMP,1
          WRITE(6,2520)XORA-40.0D0,YOR-80.0D0
          WRITE(6,2570)0.0D0
          WRITE(6,2520)XORA-40.0D0,YOR+YAMP-250.0D0
          WRITE(6,2580)ARHMAX
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,,12HP(RAD/SEC)<<)
 2560     FORMAT(3HX2,,F6.1,1H,I2,2H<<)
 2570     FORMAT(9HR900,S30,,1HP,F4.1,2H<<)
 2580     FORMAT(9HR900,S30,,1HP,D11.3,2H<<)
C
        END IF
C
C     * AMPLITUDE SPECTRA DATA PLOT ( W ) *
C
        DO 300 MODE=NSMD,NCFOM
C
          WRITE(*,2501)MODE
          READ(*,*)NLIN
 2501     FORMAT(/' Line Type for Vertical Amp.Spec.(ARV) : Mode =',I5,
     *           /'   : Solid (0) or Dot (1) ? :',\)
          WRITE(6,2502)
          IF(NLIN.EQ.1)WRITE(6,2503)
 2502     FORMAT(2HL0)
 2503     FORMAT(3HL3,,3HB30)
          WRITE(*,2504)
          READ(*,'(A)')NP
 2504     FORMAT(/' Pen Number (J1-J9) ? :',\)
          WRITE(6,2505)NP
 2505     FORMAT(A)
C
          WRITE(6,2700)XORA+CFOM(MODE)*XFAC,YOR
          DO 310 IOM=1,NOM
            OMX=DOM*DFLOAT(IOM)
            IF(OMX.LT.CFOM(MODE))GOTO 310
            WRITE(6,2710)XORA+OMX*XFAC,YOR+DABS(ARH(IOM,MODE)*AMPFAC)
            IF(OMX.GE.XVAL)GO TO 300
  310     CONTINUE
C
 2700     FORMAT(1HM,F6.1,1H,F6.1,1H,)
 2710     FORMAT(1HD,F6.1,1H,F6.1,1H,)
C
  300   CONTINUE
C
C    ** OUTPUT OF AMPLITUDE U SPECTRA PLOTTER DATA **
C
        IF(IYAMP.EQ.0)THEN
          WRITE(6,2800)
          GOTO 50
        END IF
 2800   FORMAT(2HJ2)
C
        XORA=XORA+XVAL*XFAC+200.D0
        IF((NAXY.NE.0).AND.(NVH.NE.0))THEN
C
C     * PEN, TERM, AXIS AND CHARACTER *
C
      
          WRITE(6,3000)
          WRITE(6,3010)
          WRITE(6,3020)XORA,YOR
          WRITE(6,3030)XLEN,NUX
          WRITE(6,3020)XORA-110.0D0,YOR-60.0D0
          WRITE(6,3040)0.0D0
          WRITE(6,3020)XORA+XLEN-100.0D0,YOR-60.0D0
          WRITE(6,3040)XVAL
          WRITE(6,3020)XORA+XLEN-220.0D0,YOR-120.0D0
          WRITE(6,3050)
          WRITE(6,3020)XORA,YOR
          WRITE(6,3060)YAMP,1
          WRITE(6,3020)XORA-40.0D0,YOR-80.0D0
          WRITE(6,3070)0.0D0
          WRITE(6,3020)XORA-40.0D0,YOR+YAMP-250.0D0
          WRITE(6,3080)ARHMAX
C
 3000     FORMAT(2HJ1)
 3010     FORMAT(3H=<<)
 3020     FORMAT(1HM,F6.1,1H,F6.1,1H,)
 3030     FORMAT(3HX3,,F6.1,1H,I2,2H<<)
 3040     FORMAT(7HR0,S30,,1HP,F5.1,2H<<)
 3050     FORMAT(7HR0,S30,,12HP(RAD/SEC)<<)
 3060     FORMAT(3HX2,,F6.1,1H,I2,2H<<)
 3070     FORMAT(9HR900,S30,,1HP,F4.1,2H<<)
 3080     FORMAT(9HR900,S30,,1HP,D11.3,2H<<)
C
        END IF
C
C     * AMPLITUDE SPECTRA DATA PLOT ( U ) *
C
   50   DO 500 MODE=NSMD,NCFOM
C
          WRITE(*,3001)MODE
          READ(*,*)NLIN
 3001     FORMAT(/' Line Type for Horizontal Amp.Spec.(ARH) : Mode =',I5,
     *           /'   : Solid (0) or Dot (1) ? :',\)
          WRITE(6,3002)
          IF(NLIN.EQ.1)WRITE(6,3003)
 3002     FORMAT(2HL0)
 3003     FORMAT(3HL3,,3HB30)
          WRITE(*,3004)
          READ(*,'(A)')NP
 3004     FORMAT(/' Pen Number (J1-J9) ? :',\)
          WRITE(6,3005)NP
 3005     FORMAT(A)
          WRITE(6,3210)XORA+CFOM(MODE)*XFAC,YOR
C
          DO 510 IOM=1,NOM
            OMX=DOM*DFLOAT(IOM)
            IF(OMX.LT.CFOM(MODE))GOTO 510
            ARHU=DABS(ARH(IOM,MODE)*AMPFAC*UW0(IOM,MODE))
            WRITE(6,3220)XORA+OMX*XFAC,YOR+ARHU
            IF(OMX.GE.XVAL)GO TO 500
  510     CONTINUE
C
 3210     FORMAT(1HM,F6.1,1H,F6.1,1H,)
 3220     FORMAT(1HD,F6.1,1H,F6.1,1H,)
C
  500   CONTINUE
C
C     * SKIP OF PLOTTER PEN *
C
        WRITE(6,3300)XORA+XVAL*XFAC+500.D0,YOR
 3300   FORMAT(1HM,F6.1,1H,F6.1,1H,)
C
C *** POTTER DATA FILE CLOSE ***
C
      CLOSE(6)
C
C *** AXIS DATA OUTPUT TO FILE ***
C
      YOR=YORO
      IF(NAXY.NE.0)THEN
        OPEN(5,FILE='AXH.DAT')
        WRITE(5,3500)XOR,YOR,XFAC,YFAC,AMPFAC,XVAL,IYAMP,YLEN,NVH
        CLOSE(5)
 3500   FORMAT(6D15.7,I5,D15.7,I5)
      END IF
C
      WRITE(*,3700)
 3700 FORMAT(/' USE FOR PLOTTER : A> COPY PH.DAT PRN')
      WRITE(*,*)
C
      RETURN
      END
