C     ******************************************************************
C     ******************************************************************
C     * PLOT PROGRAM FOR FOR3D MODEL.                                  *
C     * PLOTS PROPAGATION LOSS VS RANGE ON LN03 LASER PLOTTER          *
C     * SEE MAIN PROGRAM FOR3D FOR DEFINITIONS OF FOR3D VARIABLES.     *
C     ******************************************************************
C     ******************************************************************
C     * INPUT                                                          *
C     ******************************************************************
C     * INPUT UNIT NUMBER = NIU                                        *  
C     * INPUT FILE NAME   = PLT3D.IN                                   *
C     * CONTENTS:  INPUT RUNSTREAM IN FREE FORMAT                      *
C     ******************************************************************
C     * LINE 1: PHI,Z1,Z2,Z3,Y1,Y2,Y3,YL,X1,X2,X3,XL,FACT,XAVG         *
C     *   .   :   .         .         .         .         .            *
C     * LINE N:   .         .         .         .         .            *
C     ******************************************************************
C     * QUICK REFERENCE AND NOTES FOR CARD INPUT                       *
C     ******************************************************************
C               PLOT PROPAGATION LOSS VS RANGE.
C         PHI = RELATIVE BEARING OF SOLUTION - PORT SIDE IS NEGATIVE
C         Z1 = FIRST RECEIVER DEPTH TO PLOT - METERS
C         Z2 = LAST RECEIVER DEPTH TO PLOT  - METERS
C         Z3 = RECEIVER DEPTH INCREMENT - METERS
C
C         Y1 = LABEL OF Y-AXIS AT ORIGIN  - DB
C         Y2 = LABEL AT TOP OF Y-AXIS - DB
C         Y3 = INCREMENT OF Y-AXIS LABELS - DB
C         YL = LENGTH OF Y-AXIS IN INCHES
C
C         X1 = LABEL OF X-AXIS AT ORIGIN - KM - USUALLY 0
C         X2 = LABEL AT RIGHT OF X-AXIS - KM - USUALLY MAXIMUM RANGE
C         X3 = INCREMENT OF X-AXIS LABELS - KM - RANGE INCREMENT
C         XL = LENGTH OF X-AXIS IN INCHES
C
C         FACT = SCALE OF PLOT: 1.0 = FULL SIZE ; .5 = 1/2 SIZE ; ETC.
C
C         XAVG = RANGE OVER WHICH TO COMPUTE RUNNING AVERAGE -  METERS
C                IF XAVG = 0, ALL POINTS ARE PLOTTED
C     ******************************************************************
C     ******************************************************************
      PARAMETER MAXP=20000,MXPHI=360
      PARAMETER MXLYR=102,MXN=10000,NIU=1,NOU=2,NPU=6,PLTU=3
      COMPLEX HNK,HNKL,CTEMP,U(MXN),CY
      DIMENSION BETA(MXLYR,MXPHI),RHO(MXLYR,MXPHI),ZLYR(MXLYR,MXPHI),
     CIBUF(2000)
      DIMENSION P(MAXP),R(MAXP),YY(MAXP)
      DIMENSION UR(MAXP),UI(MAXP)
      DATA PI/3.141592654/,DEG/1.0/
      DATA CNVKM/1000.0/
      IPRNT=0
C
C     *** ASSIGN FOR3D OUTPUT DATA FILE
      CALL ASSIGN(NOU,'FOR3D.OUT')
C
C     *** ASSIGN PLOT INPUT PARAMETERS FILE
      CALL ASSIGN(NIU,'PLT3D.IN')
C
C     *** READ PLOT PARAMETERS
100   READ(NIU,*,END=999)PHI,Z1,Z2,Z3,Y1,Y2,Y3,YL,X1,X2,X3,XL,
     CFACT,XAVG
      IF(FACT.LE.0.0) FACT=1.0
      YINC=(Y2-Y1)/YL
      IF(Z3.EQ.0.0) Z3=1.0
C
C     *** INITIALIZE LN03 PLOTTER
      CALL PLOTS(0,0,1)
C
C     *** DEFINE ORIGIN 
      CALL PLOT (0.0,0.5,-3)
C
C     *** SET SCALE FACTOR
      CALL FACTOR(FACT)
C
C     *** GENERATE LOSS VS RANGE PLOT
      DO 270 ZR=Z1,Z2,Z3
      ZRR=ZR
      IPEN=3
      IX=1
      DX=(X2-X1)/XL
C
C     *** PLOT AXIS
      CALL AXIS2(0.,0.,'RANGE (KM)',-10,XL,0.,X1,X3,X2)
      CALL AXIS2(XL,0.,' ',-1,YL,90.,Y1,Y3,Y2)
      CALL AXIS2(0.,0.,'PROPLOSS (DB)',+13,YL,90.,Y1,Y3,Y2)
      CALL PLOT(0.0,YL,3)
      CALL PLOT(XL,YL,2)
C
C     *** READ INITIAL FOR3D PARAMETERS 
      REWIND(NOU)
      READ(2,994)NDIM,FRQ,ZS,C0,ISF,R0,Z0,N,IHNK,ITYPES,ITYPEB,ITYPPW,
     CITYPSW,FLDW,NSEC,NSOL,RMAX,DR,WDR,WZ1,WZ2,WDZ,DZ,
     CDOUGRA,NDIV,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12
994   FORMAT(1X,I1,3F10.2,I3,2F8.0,I5,5I2,F9.3,/,2I5,F10.1,F5.1,
     C5F7.1,F10.1,I4,/,(5F15.7))
C
      NPHI=NSEC
      IF(NDIM.EQ.1)THEN
        DPHI=0.0
      ELSE
        DPHI=FLDW/NSEC
      ENDIF
C
      IF(XAVG.LT.WDR) XAVG=WDR
      L=0
115   CONTINUE
      RAVG=0.0
      PLAVG=0.0
C
C     *** READ SOLUTION FIELD
      GO TO 120 ! FIRST PASS
121   WRITE(6,*)'ERROR',ANG,NN,RA,WDZ1
120   READ(NOU,995,END=170,ERR=121)ANG,NN,RA,WZ1,WDZ1,(U(I),I=1,NN)
995   FORMAT(1X,F10.5,I6,F8.0,F8.1,F8.1,/,2(E15.7,E15.7))
      IF(RA.LE.0.0) GO TO 120
      IF(ABS(ANG-PHI/DEG).GT.0.01)GO TO 120
      HNK=HNKL(2.0*PI*FRQ*RA/C0)
      I=0
      INTERP=0
      I=(ZRR-WZ1)/WDZ+1
      IF(I.GT.NN) GO TO 115
      IF(I.GE.1) GO TO 130
      I=1
      ZRR=WZ1
130   IF(WZ1+(I-1)*WDZ.NE.ZRR) INTERP=1
      IF(I.EQ.NN.AND.INTERP.EQ.1) GO TO 115
      Y=CABS(U(I))
      IF(IHNK.NE.0) Y=CABS(U(I)*HNK)
      CTEMP=U(I)+(U(I+1)-U(I))*(ZRR-(WZ1+(I-1)*WDZ))/WDZ
      IF(IHNK.NE.0) Y=CABS(CTEMP*HNK)
      IF(IHNK.EQ.0) Y=CABS(CTEMP)
      IF(Y.LE.0.0) GO TO 120
      L=L+1
      P(L)=Y
      R(L)=RA
      GO TO 120
C
170   CONTINUE
      K=0
200   K=K+1
      RAVG=0
      PLAVG=0
      NAVG=0
      DO 210 J=K,L
      IF(R(J)-R(K).GE.XAVG) GO TO 220
      RAVG=RAVG+R(J)
      PLAVG=PLAVG+P(J)
      NAVG=NAVG+1
210   CONTINUE
220   CONTINUE
      IF(NAVG.EQ.0) GO TO 250
      BIAS=0.0
      RA=RAVG/NAVG
      IF(IHNK.EQ.0.AND.RA.GT.0.0) BIAS=10.0*ALOG10(RA)
      Y=PLAVG/NAVG
      IF(Y.LE.0.0) GO TO 200
      Y=-20.0*ALOG10(Y)+BIAS
      TEMP=RA/1000.0
      IF(RA/CNVKM.GT.X2-XAVG/CNVKM) GO TO 250
      IF(RA/CNVKM.LT.X1) GO TO 200
      X=(RA/CNVKM-X1)/DX
      Y=(Y-Y1)/YINC
      IF(Y.LT.0.0) Y=0.0
      IF(Y.GT.YL) Y=YL
      CALL PLOT(X,Y,IPEN)
      IPEN=2
      GO TO 200
C
C     *** PLOT PARAMETERS
250   CONTINUE
      CALL BLOCK(XL,YL,FRQ,ZS,C0,ISF,R0,Z0,N,IHNK,ITYPEB,ITYPES,
     CRMAX,DR,WDR,WDZ,DZ,ZRR,ZLYR,BETA,RHO,XAVG,PHI,
     CDPHI,FLDW,NDIM,DOUGRA,NDIV)
      CALL PLOT(0.0,0.0,-3)
      CALL PLOT(0.0,-0.5,-999)
270   CONTINUE
      GO TO 100
C 
999   CONTINUE
      CALL PLOT(0.0,-0.5,999)
      STOP
      END
      SUBROUTINE BLOCK(XL,YL,FRQ,ZS,C0,ISF,R0,Z0,N,IHNK,ITYPEB,ITYPES,
     CRMAX,DR,WDR,WDZ,DZ,ZRR,ZLYR,BETA,RHO,XAVG,PHI,
     CDPHI,FLDW,NDIM,DOUGRA,NDIV)
      PARAMETER MXLYR=101,MXPHI=360
      DIMENSION ZLYR(MXLYR,MXPHI),BETA(MXLYR,MXPHI),RHO(MXLYR,MXPHI)
      CHARACTER DATIM(18)
      CALL DATE(DATIM)
      CALL TIME(DATIM(11))
	write(6,*)'wdz ',wdz
      DATIM(10)=' '
      NC=30 ! MAX CHAR IN STRING
      HT=.1
      DY=1.5*HT
      XBLK=5.0*HT
      YBLK=YL+(26)*DY
      CALL SYMBOL(XBLK,YBLK,HT,'FOR3D MODEL',0.,11)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'CRAY VERSION',0.,12)
      YBLK=YBLK-DY
      IF(NDIM.eq.1)CALL SYMBOL(XBLK,YBLK,HT,'2D SOLUTION',0.,11)
      IF(NDIM.EQ.2)CALL SYMBOL(XBLK,YBLK,HT,'N X 2D SOLUTION',0.,15)
      IF(NDIM.eq.3)CALL SYMBOL(XBLK,YBLK,HT,'3D SOLUTION',0.,11)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'INITIAL PARAMETERS',0.0,18)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'FRQ = ',0.,6)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,FRQ,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' HZ',0.,3)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'ZS = ',0.,5)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,ZS,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'C0 = ',0.,5)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,C0,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M/SEC',0.,6)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'R0 = ',0.,5)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,R0,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'Z0 = ',0.,5)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,Z0,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'N =  ',0.,5)
      FN=N
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,FN,0.,-1)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'DOUGRA = ',0.,9)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,DOUGRA,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'NDIV =  ',0.,8)
      FN=NDIV
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,FN,0.,-1)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'DR = ',0.,5)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,DR,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'WDR = ',0.,6)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,WDR,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'RMAX = ',0.,7)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,RMAX,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'DZ = ',0.,5)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,DZ,0.,2)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'WDZ = ',0.,6)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,WDZ,0.,2)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'DPHI = ',0.,7)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,DPHI,0.,3)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' DEGREES',0.,8)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'FLDW = ',0.,7)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,FLDW,0.,3)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' DEGREES',0.,8)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'ISF = ',0.,6)
      FPN=ISF
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,FPN,0.,-1)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'IHNK = ',0.,7)
      FPN=IHNK
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,FPN,0.,-1)
      CALL WHERE(X999,Y999,RFACT)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'ITYPES = ',0.,9)
      FPN=ITYPES
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,FPN,0.,-1)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'ITYPEB = ',0.,9)
      FPN=ITYPEB
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,FPN,0.,-1)
600   CONTINUE
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'AVG =  ',0.,7)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,XAVG,0.,-1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
601   YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'RECEIVER DEPTH = ',0.,17)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,YBLK,HT,ZRR,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,YBLK,HT,' M',0.,2)
      YBLK=YBLK-DY
      CALL SYMBOL(XBLK,YBLK,HT,'RELATIVE BEARING = ',0.,19)
      CALL WHERE(X999,Y999,RFACT)
      CALL NUMBER(X999,Y999,HT,PHI,0.,1)
      CALL WHERE(X999,Y999,RFACT)
      CALL SYMBOL(X999,Y999,HT,' DEGREES',0.,8)
      YBLK=YBLK-DY
      RETURN
      END