C
C *** READ OF ONE CHANNNEL WAVE DATA FROM DISK FILE ***
C
      SUBROUTINE WVRD(TITL,WV,NWV,MWV,DT)
C
C TITL : TITLE OF EACH CHANNEL WAVE
C WV   : WAVE AMPLITUDE DATA
C NWV  : WAVE DATA NUMBER
C MWV  : NWV IN MAIN PROGRAM DIMENSION
C DT   : DELTA TIME (SEC) 
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*30, TITL*50
      DIMENSION WV(MWV)
      INTEGER*2 IWV(128)
C
C *** INPUT OF WAVE DATA NAME AND CHANNEL NUMBER ***
C
      WRITE(*,*)
      WRITE(*,1000)
   10 WRITE(*,1010)
      READ(*,'(A)',ERR=10)FNAME
      IF (FNAME(1:1).EQ.' ') GO TO 10
   12 WRITE(*,1020)
      READ(*,*,ERR=10)ICH
 1000 FORMAT(' FROM SUBROUTIME LVRD :')
 1010 FORMAT(/' Enter wave file name :',\)
 1020 FORMAT(/' Enter output channel :',\)
C
C *** INPUT OF WAVE DATA FROM FILE DATA ***
C
      OPEN (5,FILE=FNAME,FORM='UNFORMATTED',STATUS='OLD',
     &      ACCESS='DIRECT',RECL=256,ERR=10)
C
C   * READ OF CHANNEL-POINTER ( RECORD INFORMATION ) *
C
      READ(5,REC=1) IWV
C
C   * READ OF TITL *
C
      LREC = IWV(ICH)
      READ(5,REC=LREC,ERR=20) TITL
      WRITE(*,*)TITL
      WRITE(*,1100)
      READ(*,*)NYES
      IF(NYES.NE.0)THEN
        CLOSE(5)
        GOTO 12
      END IF
 1100 FORMAT(/' Is this ch. ok ? (yes = 0) :',\)
C
C   * READ OF DATA NUMBER, DELTA TIME AND AMPLITUDE FACTER *
C
      LREC = LREC+1
      READ(5,REC=LREC) IWV
      NWV=IWV(1)
      NDT=IWV(3)
      NAF=IWV(5)
      NDTFLG=IWV(8)
      NOD=IWV(10)
      IOD=IWV(12)
      XNOD=FLOAT(NOD)
      IF(NDTFLG.NE.1) THEN
        DT = 1.0/FLOAT(NDT)
        WRITE(*,*) 'NDTFLG=',NDTFLG
      ELSE
        DT=XNOD*(10.0**IOD)
      ENDIF
      AFAC= 2.0D0**NAF
C
C   * DATA ECHO END ENTER OF DELTA TIME SKIP *
C
      NWVO=NWV
      DTO=DT
      NSKIP=1
   13 WRITE(*,1200)NWV,DT
C
   14 WRITE(*,1210)
      READ(*,*)ISYS
      IF(ISYS.EQ.0)THEN
        WRITE(*,1220)
        READ(*,*)NSKIP
        NSKIP=NSKIP+1
        NWV=NWVO/NSKIP
        DT=DTO*DFLOAT(NSKIP)
        WRITE(*,1200)NWV,DT
        WRITE(*,1230)
        READ(*,*)NSY
        IF(NSY.NE.0)GOTO 14
      END IF
 1200 FORMAT(/' Total wave data num. =',I7,' ,Delta time (sec) =',F5.3)
 1210 FORMAT( ' Skip of delta time ? (yes = 0) : ',\)
 1220 FORMAT( ' Enter skip number : ',\)
 1230 FORMAT(/' Are these skip data OK ? (yes=0) : ',\)
C
C   * INPUT OF WAVE DURATION *
C
      DUR=DFLOAT(NWV-1)*DT
      WRITE(*,1300)DUR
      WRITE(*,1310)
      READ(*,*)IYES
      IF(IYES.EQ.0)GOTO 16
      WRITE(*,1320)
      READ(*,*)DUR
      NWV=IDINT(DUR/DT)+1
      WRITE(*,1330)NWV
 1300 FORMAT(/' Total duration (sec) =',F10.5)
 1310 FORMAT(' Total wave output ? (yes = 0) :',\)
 1320 FORMAT(' Enter output duration (sec) :',\)
 1330 FORMAT(' New wave data num. =',I7)
C
C *** DATA NUMBER CHECK ***
C
   16 IF(NWV.GT.MWV)THEN
        WRITE(*,*)
        WRITE(*,*)' NWV > MWV, thus you need data skip !!!'
        WRITE(*,*)' NWV, MWV =',NWV,MWV
        WRITE(*,*)
        GOTO 13
      END IF
C
C   * READ OF WAVE DATA *
C
      NN = 128
      NREC = NWVO/NN+1
      LREC = LREC+1
      NC=0
      DO 110 J = 1,NREC
        READ(5,REC=LREC) IWV
        NUM0= NN*(J-1)
        NUM1= NUM0+NSKIP-1
        DO 100 I = 1,NN
          IF(NUM0+I.GT.NWVO)GOTO 120
          IF(MOD(NUM1+I,NSKIP).NE.0)GOTO 100
          NC=NC+1
          WV(NC) = DFLOAT(IWV(I))*AFAC
  100   CONTINUE
        LREC = LREC+1
  110 CONTINUE
C
  120 CLOSE(5)
C
      WRITE(*,1400)
      READ(*,*)NOK
      IF(NOK.NE.0)GOTO 10
      WRITE(*,*)
 1400 FORMAT(/' All data OK ? (yes=0) : ',\)
C
      RETURN
C
C *** ERROR ***
C
   20 WRITE(*,*) 'Channel access error.'
      STOP
C
      END
C
C
C
C *** FOURIER TRANSFORM FROM WAVE TO FOURIER COEFFICIENTS
C *** USE SUBROUTINE WVFAST
C
      SUBROUTINE WVFR(WV,NWV,MWV,FOUR,NFAST,DT)
C
C INPUT DATA
C WV(I)   : WAVE DATA
C NWV     : WAVE DATA NUMBER
C MWV     : NWV IN MAIN PROGRAM DIMENDION
C OUTPUT DATA
C FOUR(I) : FOUREIR COEFFICIENTS
C          (FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0)
C NFAST   : DATA NUMBER OF FOURIER SPECTRUM (2**N) > NWV
C DT      : DELTA TIME (SEC)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION WV(MWV)
      COMPLEX*16 FOUR(MWV)
C
C *** SEARCH OF FOURIER SPECTRUM DATA NUMBER ( = 2**IFAST ) > NWN ***
C
      NFAST=2
      DO 100 I=2,15
        NFAST=2*NFAST
        IF(NFAST.GE.NWV)GOTO 10
  100 CONTINUE
      WRITE(*,*)
      WRITE(*,*)' FROM SUBROUTINE WVFR :'
      WRITE(*,*)' WAVE DATA NUMBER > 2**15(=32768), THUS STOP !!!'
      STOP
   10 IFAST=I
C
C *** DATA ECHO AND DURATION ***
C
      WRITE(*,*)
      WRITE(*,*)'FROM SUBROUTINE WVFR :'
   20 WRITE(*,1000)NFAST,IFAST
      WRITE(*,1010)DFLOAT(NFAST)*DT
      WRITE(*,1020)
      READ(*,*)IDUR
 1000 FORMAT( ' TOTAL WAVE DATA NUMBER IN FFT (2**N) =',I6,',  N =',I3)
 1010 FORMAT( ' DURATION (SEC)  =',F10.5)
 1020 FORMAT(/' DOUBLE DURATION (YES=0) ?:',\)
      WRITE(*,*)
C
      IF(IDUR.EQ.0)THEN
        IFAST=IFAST+1
        NFAST=NFAST*2
        GOTO 20
      END IF
C
C *** MWV CHECK ***
C
      IF(NFAST.GT.MWV)THEN
        WRITE(*,*)
        WRITE(*,*)' FROM SUBROUTINE WVFR :'
        WRITE(*,*)' NFAST > MWV, THUS STOP !!!'
        WRITE(*,*)' NFAST, MWV =',NFAST,MWV
        STOP
      ENDIF
C
C *** INITIALIZE OF WAVE DATA
C
      DO 200 I=1,NWV
        FOUR(I)=DCMPLX(WV(I),0.0D0)
  200 CONTINUE
      DO 300 I=NWV+1,NFAST
        FOUR(I)=(0.0D0,0.0D0)
  300 CONTINUE
C
C *** FOURIER TRANSFORM FROM WAVE TO SPECTRUM ***
C
      CALL WVFT(NFAST,FOUR,MWV,-1)
C
      DFAST=DFLOAT(NFAST)
      DO 400 I=1,NFAST
        FOUR(I)=FOUR(I)/DFAST
  400 CONTINUE
C
      RETURN
      END
C
C
C
C
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE WVFT(N,A,MOM,IND)
C
C N   : NUMBER OF FOURIER COEFF. OR WAVE ( = 2**I )
C A(I): FOURIER COEFFICIENTS ( N/2+1 ) OR WAVE ( REAL PART )
C     ( FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0 )
C MOM : N IN MAIN DIMMENSION
C IND : = -1 : TRANS.,   = 1 : INVERSE TRANS.
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 A(MOM),TEMP,THETA
C
C *** CALUCURATION OF FFT ***
C
      PAI=3.1415926540D0
      J=1
      DO 140 I=1,N
        IF(I.GE.J) GO TO 110
        TEMP=A(J)
        A(J)=A(I)
        A(I)=TEMP
  110   M=N/2
  120   IF(J.LE.M) GO TO 130
        J=J-M
        M=M/2
        IF(M.GE.2) GO TO 120
  130   J=J+M
  140 CONTINUE
      KMAX=1
  150 IF(KMAX.GE.N) RETURN
      ISTEP=KMAX*2
      DO 170 K=1,KMAX
        THETA=DCMPLX(0.0D0,PAI*DFLOAT(IND*(K-1))/DFLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*CDEXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END
C
C
C
C
C *** OUTPUT OF WAVE DATA TO FILE ***
C
      SUBROUTINE WVOT(WV,NWV,MWV,DT)
C
C WV(I) : WAVE DATA
C NWV   : NUMBER OF WAVE DATA
C MWV   : NWV IN MAIN
C DT    : DELTA TIME
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION WV(MWV)
      CHARACTER FNAME*30
C
C *** FILE OUTPUT TO FILE DATA ***
C
      WRITE(*,1000)
      WRITE(*,1010)
      READ(*,*)IYES
      IF(IYES.NE.0)GOTO 20
 1000 FORMAT(/' From SUBROUTINE WVOT : wave data to file')
 1010 FORMAT(/' File out (yes=0)? :',\)
C
   10 WRITE(*,1020)
      READ(*,'(A)',ERR=10)FNAME
      IF(FNAME(1:1).EQ.' ')GOTO 10
 1020 FORMAT(/' Enter file name :',\)
C
      OPEN(7,FILE=FNAME,ERR=10)
      WRITE(7,2000)NWV,DT
      DO 100 IWV=1,NWV,8
      WRITE(7,2010)(WV((IWV-1)+I),I=1,8)
  100 CONTINUE
      CLOSE(7)
 2000 FORMAT(I5,3F10.4)
 2010 FORMAT(8F10.4)
C
   20 WRITE(*,*)
C
      RETURN
      END
C
C
C
C
C *** OUTPUT OF AMPLITUDE SPETRUM TO PLOTTER ***
C
      SUBROUTINE SPPT(FOM,DOM,NOM,MOM,SPCA)
C
C FOM     : FIRST CIRCULAR FREQUENCY : OMEGA OF SPCA(1)
C DOM     : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C NOM     : TOTAL NUMBER OF CIRCULAR FREQUENCY
C MOM     : NOM IN MAIN PROGRAM DIMENSION
C SPCA(I) : FOURIER SPECTRUM AMPLITUDE
C  ( FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0 )
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION SPCA(MOM)
      CHARACTER YUNIT*10
C
C *** SEARCH OF MAXIMUM AMPLITUDE AND FINAL OMEGA ***
C
      AMX=0.0D0
      DO 100 IOM=1,NOM
        APC=SPCA(IOM)
        IF(APC.GE.AMX)THEN
        AMX=APC
        IOMMX=IOM
        END IF
  100 CONTINUE
C
      EOM=FOM+DOM*DFLOAT(NOM-1)
      OMMX=FOM+DOM*DFLOAT(IOMMX-1)
C
C *** OUTPUT OF DATA INFORMATION ***
C
      WRITE(*,*)
      WRITE(*,1000)
   10 WRITE(*,1010)
      WRITE(*,1020)DOM,NOM
      WRITE(*,1030)FOM,EOM
      WRITE(*,1040)AMX,OMMX
 1000 FORMAT(' FROM SUBROUTINE SPPT : AMPLITUDE SPETRUM PLOT')
 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(' MAXIMUM AMPLITUDE = ',F10.5,' , OMEGA = ',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='AX.DAT')
        READ(5,1210)XOR,YOR,XFAC,YFAC,OLDFOM,XVAL,NST0
        CLOSE(5)
 1210   FORMAT(6D15.7,I5)
        OM1=FOM
        IF(NST0.NE.0)OM1=FOM-OLDFOM
C
      ELSE
C
        WRITE(*,1300)
        WRITE(*,1310)
        READ(*,*)NST0
        OM1=0.0D0
        OM2=FOM
        IF(NST0.EQ.0)OM1=FOM
        IF(NST0.EQ.0)OM2=0.0D0
        WRITE(*,1320)
        READ(*,*)XVAL
        WRITE(*,1330)
        READ(*,*)XLEN
        XLEN=XLEN*100.0D0
        WRITE(*,1340)
        READ(*,*)NUX
        XFAC=XLEN/(XVAL-OM2)
        WRITE(*,1350)
        READ(*,*)YVAL
        WRITE(*,1360)
        READ(*,*)YLEN
        YLEN=YLEN*100.0D0
        WRITE(*,1370)
        READ(*,*)NUY
        YFAC=YLEN/YVAL
        WRITE(*,1380)
        READ(*,*)IYUT
        IF(IYUT.EQ.1)YUNIT='  (CM*SEC)'
        IF(IYUT.EQ.2)YUNIT='(KINE*SEC)'
        IF(IYUT.EQ.3)YUNIT=' (GAL*SEC)'
        IF(IYUT.EQ.4)YUNIT='          '
 1300   FORMAT(' AXIS DATA INPUT :')
 1310   FORMAT(' X-AXIS START FROM 0.0 OR FIRST OMEGA ?'
     *        /' IF START 0.0 OMEGA, PLEASE ENTER 0 :',\) 
 1320   FORMAT(/' ENTER X-AXIS MAX (MAX OMEGA)                  :',\)
 1330   FORMAT( ' ENTER X-AXIS LENGTH (OMEGA: cm)               :',\)
 1340   FORMAT( ' ENTER X-AXIS UNIT NUMBER                      :',\)
 1350   FORMAT(/' ENTER Y-AXIS MAX VALUE (Y-AXIS AMP > MAX AMP) :',\)
 1360   FORMAT( ' ENTER Y-AXIS LENGTH (AMP : cm)                :',\)
 1370   FORMAT( ' ENTER Y-AXIS UNIT NUMBER                      :',\)
 1380   FORMAT(/' 1: DISP(CM*SEC),  2: VEL(KINE*SEC),'/,
     *         ' 3: ACC(GAL*SEC),  4: NON DIMENDION :',\)
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 = 0) :',\)
C
      END IF
C
C *** OUTPUT PLOTTER DATA ***
C 
      OPEN(6,FILE='PS.DAT')
C
      IF(NAXY.EQ.0)GOTO 20
C
C   * PEN, TERM, AXIS AND CHARACTER *
C
      WRITE(6,2000)
      WRITE(6,2010)
      WRITE(6,2020)XOR,YOR
      WRITE(6,2030)XLEN,NUX
      WRITE(6,2020)XOR-110.D0,YOR-60.D0
      WRITE(6,2040)OM2
      WRITE(6,2020)XOR+XLEN-100.D0,YOR-60.D0
      WRITE(6,2040)XVAL
      WRITE(6,2020)XOR+XLEN-220.D0,YOR-120.D0
      WRITE(6,2050)
      WRITE(6,2020)XOR,YOR
      WRITE(6,2060)YLEN,NUY
      WRITE(6,2020)XOR-40.D0,YOR-80.D0
      WRITE(6,2070)0.D0
      WRITE(6,2020)XOR-40.D0,YOR+YLEN-80.D0
      WRITE(6,2075)YVAL
      WRITE(6,2020)XOR-100.D0,YOR+YLEN-250.D0
      WRITE(6,2080)YUNIT
 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<<)
 2075 FORMAT(9HR900,S30,,1HP,F6.1,2H<<)
 2080 FORMAT(9HR900,S30,,1HP,A10,2H<<)
C
C   * DATA PLOT *
C
   20 WRITE(6,2100)XOR+OM1*XFAC,YOR+SPCA(1)*YFAC
      DO 200 IOM=2,NOM
      OMX=OM1+DOM*DFLOAT(IOM-1)
      WRITE(6,2110)XOR+OMX*XFAC,YOR+SPCA(IOM)*YFAC
      IF(OMX.GE.XVAL)GO TO 30
  200 CONTINUE
   30 WRITE(6,2100)XOR+(OM1+DOM*DFLOAT(IOM))*XFAC+500.D0,YOR
 2100 FORMAT(1HM,F6.1,1H,F6.1,1H,)
 2110 FORMAT(1HD,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='AX.DAT')
        WRITE(5,2200)XOR,YOR,XFAC,YFAC,FOM,XVAL,NST0
        CLOSE(5)
 2200   FORMAT(6D15.7,I5)
      END IF
C
      WRITE(*,3000)
 3000 FORMAT(/' USE FOR PLOTTER : A> COPY PS.DAT PRN')
      WRITE(*,*)
C
      RETURN
      END
C
C
C
C
C *** PLOT OF ONE CHANNNEL WAVE DATA FROM DISK FILE ***
C
      SUBROUTINE WVPT(TITL,WV,NWV,MWV,DT)
C
C TITL : TITLE OF EACH CHANNEL
C WV   : WAVE AMPLITUDE DATA
C NWV  : WAVE DATA NUMBER
C MWV  : NWV IN MAIN PROGRAM DIMENSION
C DT   : DELTA TIME ( SEC ) 
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER TITL*50,YUNIT*6
      DIMENSION WV(MWV)
C
C *** SEARCH OF MAXIMUM AMPLITUDE AND DURATION  ***
C
      AMX=0.0D0
      DO 100 IWV=1,NWV
      AWV=DABS(WV(IWV))
      IF(AWV.GE.AMX)THEN
      AMX=AWV
      IWVMX=IWV
      END IF
  100 CONTINUE
C
      AMX=WV(IWVMX)
      AMXTM=DT*DFLOAT(IWVMX-1)
      DUR=DT*DFLOAT(NWV-1)
C
C *** OUTPUT OF DATA INFORMATION ***
C
      WRITE(*,*)
      WRITE(*,1000)
   10 WRITE(*,1010)
      WRITE(*,1020)DUR
      WRITE(*,1030)AMX,AMXTM
 1000 FORMAT(' FROM SUBROUTINE WVPT :')
 1010 FORMAT(' DATA INFORMATION :')
 1020 FORMAT(' WAVE DURATION (SEC) = ',F8.3)
 1030 FORMAT(' MAXIMUM AMPLITUDE = ',F10.5,' ,TIME (SEC) = ',F8.3/)
C
C *** INPUT OF PLOTTER AXIS DATA ***
C
      WRITE(*,1100)
      WRITE(*,1110)
      READ(*,*)XVAL
      WRITE(*,1120)
      READ(*,*)XLEN
      XLEN=XLEN*100.0D0
      WRITE(*,1130)
      READ(*,*)NUX
      XFAC=XLEN/XVAL
      WRITE(*,1140)
      READ(*,*)YVAL
      WRITE(*,1150)
      READ(*,*)YLEN
      YLEN=YLEN*100.0D0
      WRITE(*,1160)
      READ(*,*)NUY
      YFAC=YLEN/YVAL
      WRITE(*,1170)
      READ(*,*)IYUT
      IF(IYUT.EQ.1)YUNIT='  (CM)'
      IF(IYUT.EQ.2)YUNIT='(KINE)'
      IF(IYUT.EQ.3)YUNIT=' (GAL)'
 1100 FORMAT(' CONTROLL DATA INPUT :')
 1110 FORMAT(/' ENTER X-AXIS MAX VALUE (DURATION : sec)       :',\)
 1120 FORMAT( ' ENTER X-AXIS LENGTH ON PLOTTER (cm)           :',\)
 1130 FORMAT( ' ENTER X-AXIS UNIT NUMBER                      :',\)
 1140 FORMAT(/' ENTER Y-AXIS MAX VALUE (Y-AXIS AMP > MAX AMP) :',\)
 1150 FORMAT( ' ENTER Y-AXIS LENGTH ON PLOTTER (AMP: cm)      :',\)
 1160 FORMAT( ' ENTER Y-AXIS UNIT NUMBER                      :',\)
 1170 FORMAT(/' 1:DISP(CM),2:VEL(KINE),3:ACC(GAL) ?           :',\)
C
      WRITE(*,1200)
      READ(*,*)XOR
      XOR=(XOR+2.0)*100.0D0
      WRITE(*,1210)
      READ(*,*)YOR
      YOR=(YOR+2.0D0)*100.0D0+YLEN
      WRITE(*,1220)
      READ(*,*)IYES
      IF(IYES.NE.0)GOTO 10
 1200 FORMAT(/' ENTER X LEFT SIDE EDGE POSITION (cm):',\)
 1210 FORMAT( ' ENTER Y DOWN SIDE EDGE POSITION (cm):',\)
 1220 FORMAT(/' ALL INPUT DATA OK ? (YES = 0) :',\)
C
C *** OUTPUT PLOTTER DATA ***
C 
      OPEN(6,FILE='PW.DAT')
C
C   * PEN, TERM, AXIS, UNIT CHARACTER AND TITL PLOT *
C
      WRITE(6,2000)
      WRITE(6,2010)
      WRITE(6,2020)XOR,YOR
      WRITE(6,2030)XLEN,NUX
      WRITE(6,2020)XOR+XLEN-100.0D0,YOR-60.0D0
      WRITE(6,2040)XVAL
      WRITE(6,2020)XOR+XLEN-100.0D0,YOR-120.0D0
      WRITE(6,2050)
      WRITE(6,2020)XOR,YOR-YLEN
      WRITE(6,2060)2.0D0*YLEN,2*NUY
      WRITE(6,2020)XOR-40.0D0,YOR-YLEN-60.0D0
      WRITE(6,2070)-YVAL
      WRITE(6,2020)XOR-40.0D0,YOR-80.0D0
      WRITE(6,2080)0.0D0
      WRITE(6,2020)XOR-40.0D0,YOR+YLEN-80.0D0
      WRITE(6,2070)YVAL
      WRITE(6,2020)XOR-100.0D0,YOR+YLEN-180.0D0
      WRITE(6,2090)YUNIT
      WRITE(*,2095)
      READ(*,*)ITITL
      IF(ITITL.NE.0)THEN
        WRITE(6,2020)XOR-100.0D0,YOR-YLEN-150.0D0
        WRITE(6,2100)TITL
      END IF
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,,8HP(SEC)<<)
 2060 FORMAT(3HX2,,F6.1,1H,I2,2H<<)
 2070 FORMAT(9HR900,S30,,1HP,F6.1,2H<<)
 2080 FORMAT(9HR900,S30,,1HP,F4.1,2H<<)
 2090 FORMAT(9HR900,S30,,1HP,A6,2H<<)
 2095 FORMAT(/' TITL OFF (YES=0)?:',\)
 2100 FORMAT(7HR0,S30,,1HP,A50,2H<<)
C
C   * DATA PLOT *
C
      WRITE(6,2200)XOR,YOR+WV(1)*YFAC
      DO 200 IWV=1,NWV-1
      XX=DT*DFLOAT(IWV)*XFAC
      WRITE(6,2210)XOR+XX,YOR+WV(IWV+1)*YFAC
      IF(XX.GT.XLEN)GOTO 20 
  200 CONTINUE
   20 WRITE(6,2200)XOR+DT*DFLOAT(IWV)*XFAC+500.0,YOR
 2200 FORMAT(1HM,F6.1,1H,F6.1,1H,)
 2210 FORMAT(1HD,F6.1,1H,F6.1,1H,)
C
      CLOSE(6)
C
      WRITE(*,3000)
 3000 FORMAT(' USE FOR PLOTTER : A> COPY PW.DAT PRN')
      WRITE(*,*)
C
      RETURN
      END
      