C file S3.FOR
C ----------------------------------------------------------------------!-------
! 11-17 Nov 2021 - Added APLT=5 for polar plot horizontal projection
! 14 Dec 2016 - Added labels to raypath plots using PLTNOTE
! 12 Apr 2011 - Added SUBROUTINE RAYPLUS
C 16 Mar 2010 - Increased dimension of UGP and DUGP arrays to 611
C 29 Jun 2009 - Added arguments to CIRCLE and PLOTNOTE
C 29 Jun 2009 - Split SUBROUTINE PLOTNOTE off from SUBROUTINE CIRCLE
C 29 Jun 2009 - Moved CALL CIRCLE from S0
C 25 Jun 2009 - Minor change
C 23 Jun 2009 - Changed the default character size for DDTEXT calls to 1.
C 22 Jun 2009 - Added size and bold to optional labeling
C 18 Jun 2009 - Added optional labeling of raypath plots for vortex
C 16 Jun 2009 - Added subroutine circle to draw a circle at wall of vortex
C 30 Apr 2009 - replaced REAL*8 with DOUBLE PRECISION
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 2 Mar 2006 - Moved PUTDES, NUMSTG, SFILL, STRIM, RERR,
C RERROR, STOPIT, PUTKST, PUTKCT, OPNREP, ZAPFIL, OVERRD, SFILTR, ALCOSH,
C GAUSEL, FTANH, FTANH2 to S5.for
C ----------------------------------------------------------------------!-------
Contains RAYPLT with ENTRY ENDPLT and SUBROUTINES CIRCLE, PLOTNOTE, and RAYPLUS
C ----------------------------------------------------------------------!-------
C 3 January 2006. Added SUBROUTINE FTANH2 to support 2nd derivatives.
C 6 September 2005. Added SUBROUTINE BETTERCEA.
C 17 Sept. 2004.Changed tanh to dtanh in SUBROUTINE FTANH
C 6 July 2004.  Changed May 1991 correction.
C 5 June 2004.  Separated SUBROUTINE OCNHD into its own file.
C 2 June 2004.  Added printing of pressure and molecular weight to 
C               ENTRY HEADER1 in SUBROUTINE OCNHD
C 28 April 2004. Modified SUBROUTINE FTANH.
C 25 April 2004. Modified SUBROUTINE FTANH.
C 21 April 2004. Added SUBROUTINE FTANH at end of file.
C Sometime in 2003, changed line 17 of FUNCTION RERR (ERT = ...).
C 6 May 2003. Added 3 lines at line 58 in ENTRY HEADER1.
C ----------------------------------------------------------------------!-------
      SUBROUTINE RAYPLT
C ----------------------------------------------------------------------!-------
C     MAIN PLOTTING PROGRAM; INITIALIZES, READS INPUT, PLOTS
C     PROJECTIONS OF RAYS ON A VERTICAL OR HORIZONTAL PLANE.
C        ABS(PLT)=1. PLOTS PROJECTION OF RAYPATH ON VERTICAL PLANE
C                    RECTANGULAR EXPANSION BY FACTOR 'PFACTR'
C                =2. PLOTS PROJECTION OF RAYPATH ON GROUND
C                =3. VERICAL PROJECTION USING RADIAL EXPANSION BY FACTOR
C                    'PFACTR'
C     COMMON DECK "FILEC" INSERTED HERE
 
      implicit double precision (a-h,o-z)
      CHARACTER ID*80
      COMMON /FILEC/NPLTDP
      COMMON /PLT/ XL,XR,YB,YT,PRESET
      COMMON/PLT/RMIN,RMAX,ALPHA,APLT
C
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10

C     COMMON DECK "FLAG" INSERTED HERE
      LOGICAL NEWWR,NEWWP,NEWTRC,PENET
      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH
      COMMON/FLGP/NSET

C     COMMON DECK "RKAM" INSERTED HERE
      DOUBLE PRECISION KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)
      COMMON/WW/DUM(10),MAXW,W(1000)
      common/ident/id
      DOUBLE PRECISION MAXSTP,MAXERR,INTYP,LLAT,LLON
      EQUIVALENCE (EARTHR,W(1)),(RAY,W(2)),(XMTRH,W(3)),(TLAT,W(4)),
     1 (TLON,W(5)),(OW,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),
     2 (AZ1,W(10)),(AZBEG,W(11)),(AZEND,W(12)),(AZSTEP,W(13)),
     3 (BETA,W(14)),(ELBEG,W(15)),(ELEND,W(16)),(ELSTEP,W(17)),
     8   (RUNSUP,W(18)),(RCVRH,W(20)),
     4 (ONLY,W(21)),(HOP,W(22)),(MAXSTP,W(23)),(PLAT,W(24)),(PLON,W(25))
     5,(HMAX,W(26)),(RAYFNC,W(29)),(EXTINC,W(33)),
     6 (HMIN,W(27)),(RGMAX,W(28)),
     8 (INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),
     6 (STEP1,W(44)),(STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),
     7 (SKIP,W(71)),(RAYSET,W(72)),(PRTSRP,W(74)),(HITLET,W(75))
     9 ,(BINRAY,W(76)),(PAGLN,W(77)),(PLT,W(81)),(PFACTR,W(82)),
     1 (LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))
     2,(TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))
      DOUBLE PRECISION MMODEL,MFORM,MID
C
C     WIND        100-124
      EQUIVALENCE (W(100),UMODEL),(W(101),UFORM),(W(102),UID)
      EQUIVALENCE (U0,W(103)),(R0,W(104)),(LAM0,W(105)),(PH0,W(106))    
      DOUBLE PRECISION LAM0

C     COMMON DECK "B1" INSERTED HERE
      INTEGER UMX,UNTBL,UITBL,UFRMTBL,IDSU(20)
      COMMON/B1/UMX,UNTBL(10),UITBL(10),UFRMTBL(10),UGP(611)      ! 16 Mar 2010
      EQUIVALENCE (UGP,IDSU)

C
C     DELTA WIND        125-149
      EQUIVALENCE (W(125),DUMODEL),(W(126),DUFORM),(W(127),DUID)
      EQUIVALENCE (U1,W(130)),(GAMMA1,W(131)),(LAM1,W(128)),(PH1,W(129))    

C     COMMON DECK "B2" INSERTED HERE
      INTEGER DUMX,DUNTBL,DUITBL,DUFRMTB,IDSDU(20)
      COMMON/B2/DUMX,DUNTBL(10),DUITBL(10),DUFRMTB(10),DUGP(611)  ! 16 Mar 2010
      EQUIVALENCE (DUGP,IDSDU)
C
C     SOUND SPEED 150-174
      EQUIVALENCE (W(150),CMODEL),(W(151),CFORM),(W(152),CID)
      EQUIVALENCE (W(153),REFC)
C
C     DELTA SOUND SPEED 175-199
      EQUIVALENCE (W(175),DCMODEL),(W(176),DCFORM),(W(177),DCID)
C
C     TEMPERATURE    200-224
      EQUIVALENCE (W(200),TMODEL),(W(201),TFORM),(W(202),TID)
C
C     DELTA TEMPERATURE    225-249
      EQUIVALENCE (W(225),DTMODEL),(W(226),DTFORM),(W(227),DTID)
C
C     MOLECULAR   250-274
      EQUIVALENCE (W(250),MMODEL),(W(251),MFORM),(W(252),MID)
C
C     RECEIVER HEIGHT   275-299
      EQUIVALENCE (W(275),RMODEL),(W(276),RFORM),(W(277),RID)
C
C     TOPOGRAPHY   300-324
      EQUIVALENCE (W(300),GMODEL),(W(301),GFORM),(W(302),GID)
C
C     DELTA TOPOGRAPHY   325-349
      EQUIVALENCE (W(325),GUMODEL),(W(326),GUFORM),(W(327),GUID)
C
C     ATMOSPHERIC SURFACE TOPOGRAPHY   350-374
      EQUIVALENCE (W(350),SMODEL),(W(351),SFORM),(W(352),SID)
C
C     ATMOSPHERIC SURFACE TOPOGRAPHY   375-399
      EQUIVALENCE (W(375),SUMODEL),(W(376),SUFORM),(W(377),SUID)
C     PLOT ENHANCEMENTS CONTROL PARAMETERS
C
      EQUIVALENCE (W(490),XFQMDL),(W(491),YFQMDL)
C     ABSORPTION        500-524
      EQUIVALENCE (W(500),AMODEL),(W(501),AFORM),(W(502),AID)
C
C     DELTA ABSORPTION        525-549
      EQUIVALENCE (W(525),DAMODEL),(W(526),DAFORM),(W(527),DAID)
C
C     PRESSURE        550-574
      EQUIVALENCE (W(550),PMODEL),(W(551),PFORM),(W(552),PID)
C
C     DELTA PRESSURE        575-599
      EQUIVALENCE (W(575),DPMODEL),(W(576),DPFORM),(W(577),DPID)
C
      DOUBLE PRECISION LTIC
C%      INTEGER IWDMP(411)
C%      EQUIVALENCE(IWDMP,ID)
C
      EXTERNAL PLOT,PLTLB,PLTLBN,PLTHLB
C
      DATA NWDSRK/4/
      DATA K/0/
C
      IF (.NOT.NEWWR) GO TO 8
      APLT=dABS(PLT)
C
C          NEW W ARRAY -- REINITIALIZE
      NEWWR=.FALSE.
      PRESET=1.d0
C     IF NO ACTIVE PLOTTING,WE ARE STILL DUMPING DATA
      IF(PLT.EQ.0.d0) GO TO 5
C     INITIALIZE ANNOTATION MODEL
C  CONVERT COORDINATES OF VERTICAL PLANE FROM GEOGRAPHIC TO GEOMAGNETIC
      SW=dSIN (PLAT)
      CW=dSIN (PID2-PLAT)
      SLM=dSIN (LLAT)
      CLM=dSIN (PID2-LLAT)
      SRM=dSIN (RLAT)
      CRM=dSIN (PID2-RLAT)
      CDPHI=dSIN (PID2-(LLON-PLON))
      PHL=dATAN2(dSIN (LLON-PLON)*CLM,CDPHI*SW*CLM-CW*SLM)
      CTHL=CDPHI*CW*CLM+SW*SLM
      THL=DACOS(CTHL)          ! 6 Sept 2005
      STHL=dSIN (dACOS (CTHL))
      CDPHI=dSIN (PID2-(RLON-PLON))
      PHR=dATAN2(dSIN (RLON-PLON)*CRM,CDPHI*SW*CRM-CW*SRM)
      CTHR=CDPHI*CW*CRM+SW*SRM
      THR=DACOS(CTHR)          ! 6 Sept 2005
      STHR=dSIN (dACOS (CTHR))
C      CLR=CTHL*CTHR+STHL*STHR*dSIN (PID2-(PHL-PHR))
C      SLR=DSQRT (1.0D0-CLR**2)
      CALL BETTERCEA(THL,THR,THLR,CTHL,CTHR,CLR,STHL,STHR,SLR,PHL-PHR)
C      ALPHA=.5D0*DACOS (CLR)
      ALPHA=.5D0*THLR
C
      IF (APLT.EQ.2.d0) GO TO 3
      IF (APLT.EQ.5.d0) GO TO 3                                    ! 15 Nov 2021
C
C     VERTICAL PROJECTIONS ONLY
      T=HB
      HB=DMIN1(T,HT)
      HT=DMAX1(T,HT)
      RMIN=EARTHR+HB
      YT1=YT
      RMAX=EARTHR+HT
C
      IF(APLT.NE.4.d0) GO TO 100
C
C     SCALE FOR CARTESIAN PROJECTION(=4)
      XR=ALPHA
      XR1=XR
      XL=-XR
      XL1=XL
      YB=RMIN
      YT=RMAX
      GO TO 5
C
C     PROJECTIONS 1 AND 3
100   XR1=RMIN*dSIN (ALPHA)
      XL1=-XR1
      YB=RMIN*dSIN (PID2-ALPHA)
ccccc      IF(APLT.EQ.3.d0) RMAX=RMIN+(RMAX-RMIN)*PFACTR
C.....Correction M. Jones, March,1991
ccccc      IF(APLT.EQ.3.) RMAX=earthr+(RMAX-earthr)*PFACTR
C.....new correction, May 1991
C
      if (aplt.eq.3.) rmax = rmin+(rmax-rmin)*pfactr
C     &                *(rmin/(earthr+(rmin-earthr)*pfactr))   ! 6 JULY 2004

      XR=DMAX1(RMAX*DSIN(ALPHA),(RMAX-YB)/2.D0)
      XL=-XR
      YT1=RMAX*(YB/RMIN)
      YT=2.d0*XR
      IF(APLT.EQ.1.d0) YT=YT/PFACTR
      YT=YT+YB
      GO TO 5
C
C     HORIZONTAL PROJECTION(=2)
C     PLUS POLAR PLOT HORIZONTAL PROJECTION(=5)                    ! 15 Nov 2021
3     ALPH1=dATAN2(STHR*dSIN (PHR-PHL),(CTHR-CTHL*CLR)/STHL)
      XL=0.d0
      XL1=0.d0
      XR=EARTHR*2.d0*ALPHA
      XR1=XR
C USE 90% OF X-RANGE FOR Y-RANGE
      RMAX=0.5d0*(0.9d0*XR)/PFACTR
      YT=RMAX
      RMIN=-RMAX
      YB=RMIN
      X0ANGLE=PI-THL                                               ! 15 Nov 2021
      X0PLOT=EARTHR*X0ANGLE                                        ! 15 Nov 2021
      Y0ANGLE=0.D0                                                 ! 15 Nov 2021
      Y0PLOT=EARTHR*Y0ANGLE                                        ! 15 Nov 2021
      GAM1=PID2                                                    ! 15 Nov 2021
      BET1=PID2                                                    ! 15 Nov 2021
      COSBET1=0.D0                                                 ! 15 Nov 2021
      SINBET1=1.D0                                                 ! 15 Nov 2021
      SALPH1=DSIN(ALPH1)                                           ! 15 Nov 2021
      CALPH1=DCOS(ALPH1)                                           ! 15 Nov 2021
      IF(dabs(SALPH1).lt.1.0d-5) GO TO 5                           ! 16 Nov 2021
      Y0ANGLE=-DASIN(STHL*SALPH1)                                  ! 17 Nov 2021
      Y0PLOT=EARTHR*Y0ANGLE                                        ! 15 Nov 2021
      X0ANGLE=DASIN(DTAN(Y0ANGLE)*CALPH1/SALPH1)                   ! 17 Nov 2021
      X0PLOT=EARTHR*X0ANGLE                                        ! 15 Nov 2021
      GAM1=DASIN(DTAN(Y0ANGLE)*CALPH1/(STHL*SALPH1))               ! 17 Nov 2021
      BET1=PHL-PHR-GAM1                                            ! 17 Nov 2021
      COSBET1=DCOS(BET1)                                           ! 15 Nov 2021
      SINBET1=DSIN(BET1)                                           ! 15 Nov 2021
C
5     IF(NPLTDP.LE.0) GO TO 8
      WRITE(NPLTDP) 1,NWDSRK,1,411
      WRITE(NPLTDP) id
C
8     NEW=0
      IF (NEWTRC) NEW=1
      NEWTRC=.FALSE.
C
      IF(NPLTDP.LE.0) GO TO 88
      IF(NEW.NE.1) GO TO 84
      WRITE(NPLTDP) 1,NWDSRK,17,25
      WRITE(NPLTDP) id
C
84    WRITE(NPLTDP) 3-NEW,R,TH,PH
C
88    IF(PLT.EQ.0.d0) RETURN
C
      STH=dSIN (TH)
      CTH=dSIN (PID2-TH)
C     CR=CTHR*CTH+STHR*STH*dSIN (PID2-(PHR-PH))
C     CR=1-(PHR-PH)**2/2
      CALL BETTERCEA(THR,TH,TR,CTHR,CTH,CR,STHR,STH,SR,PHR-PH)
C     CL=CTHL*CTH+STHL*STH*dSIN (PID2-(PHL-PH))
C     CL=1-(PHL-PH)**2/2
      CALL BETTERCEA(THL,TH,TL,CTHL,CTH,CL,STHL,STH,SL,PHL-PH)
      CEA=dATAN2(CR-CL*CLR,CL*SLR)
      NTEST=4.5D0+DSIGN(2.D0,CR)+DSIGN(1.D0,CL)+DSIGN(.5D0,CLR)
      SELECT CASE (NTEST)
      CASE (1)
         PRINT*,'CEA BEFORE = ',CEA
         CEA=DATAN2((-2.D0+SR**2/(1.D0-CR)+SL**2/(1.D0-CL)
     1    +SLR**2/(1.D0-CLR)-(SL**2*SLR**2/((1.D0-CL)*(1.D0-CLR))))
     2    ,CL*SLR)
         PRINT*,'CASE #: ',NTEST
         PRINT*,'CEA AFTER = ',CEA
      CASE (2)
         PRINT*,'CEA BEFORE = ',CEA
         CEA=DATAN2((SR**2/(1.D0-CR)+SL**2/(1.D0+CL)
     1    -SLR**2/(1.D0-CLR)+(SL**2*SLR**2/((1.D0+CL)*(1.D0-CLR))))
     2    ,CL*SLR)
         PRINT*,'CASE #: ',NTEST
         PRINT*,'CEA AFTER = ',CEA
      CASE (3)
         PRINT*,'CEA BEFORE = ',CEA
         CEA=DATAN2((SR**2/(1.D0-CR)+SL**2/(1.D0+CL)
     1    -SLR**2/(1.D0-CLR)+(SL**2*SLR**2/((1.D0+CL)*(1.D0-CLR))))
     2    ,CL*SLR)
c        K=K+1
c        IF (K.GT.100) THEN
c        K=0
         PRINT*,'CASE #: ',NTEST
         PRINT*,'CEA AFTER = ',CEA
c        ENDIF
      CASE (4)
         PRINT*,'CEA BEFORE = ',CEA
         CEA=DATAN2((-2.D0+SR**2/(1.D0-CR)+SL**2/(1.D0+CL)
     1    +SLR**2/(1.D0+CLR)-(SL**2*SLR**2/((1.D0+CL)*(1.D0+CLR))))
     2    ,CL*SLR)
         PRINT*,'CASE #: ',NTEST
         PRINT*,'CEA AFTER = ',CEA
      CASE (5)
         PRINT*,'CEA BEFORE = ',CEA
         CEA=DATAN2((-SR**2/(1.D0+CR)-SL**2/(1.D0-CL)
     1    -SLR**2/(1.D0-CLR)-(SL**2*SLR**2/((1.D0-CL)*(1.D0-CLR))))
     2    ,CL*SLR)
         PRINT*,'CASE #: ',NTEST
         PRINT*,'CEA AFTER = ',CEA
      CASE (6)
         PRINT*,'CEA BEFORE = ',CEA
         CEA=DATAN2((2.D0-SR**2/(1.D0+CR)-SL**2/(1.D0-CL)
     1    +SLR**2/(1.D0+CLR)+(SL**2*SLR**2/((1.D0-CL)*(1.D0+CLR))))
     2    ,CL*SLR)
         PRINT*,'CASE #: ',NTEST
         PRINT*,'CEA AFTER = ',CEA
      CASE (7)
         PRINT*,'CEA BEFORE = ',CEA
         CEA=DATAN2((2.D0-SR**2/(1.D0+CR)+SL**2/(1.D0+CL)
     1    -SLR**2/(1.D0-CLR)+(SL**2*SLR**2/((1.D0+CL)*(1.D0-CLR))))
     2    ,CL*SLR)
         PRINT*,'CASE #: ',NTEST
         PRINT*,'CEA AFTER = ',CEA
      CASE (8)
C         PRINT*,'CEA BEFORE = ',CEA
         CEA=DATAN2((-SR**2/(1.D0+CR)+SL**2/(1.D0+CL)
     1    +SLR**2/(1.D0+CLR)-(SL**2*SLR**2/((1.D0+CL)*(1.D0+CLR))))
     2    ,CL*SLR)
C         PRINT*,'CASE #: ',NTEST
C         PRINT*,'CEA AFTER = ',CEA
      END SELECT
987   FORMAT(A,E50.40)
C
      IF(APLT.NE.4.d0) GO TO 150
      CALL PLOT(CEA-ALPHA,R,NEW)
      RETURN
C
150   IF(APLT.EQ.2.d0) GO TO 10
      IF(APLT.EQ.5.d0) GO TO 160                                   ! 11 Nov 2021
      RX=R
ccccc      IF(APLT.EQ.3.d0) RX=RMIN+(R-RMIN)*PFACTR
C.....Correction M. Jones, March,1991
ccccc      if (aplt .eq. 3.0) rx = earthr + (r-earthr)*pfactr

C.....new correction, May 1991
C
      if (aplt.eq.3.) rx = rmin+(r-rmin)*pfactr
C     &               *(rmin/(earthr+(rmin-earthr)*pfactr))   ! 6 JULY 2004

      CALL PLOT(CEA-ALPHA,RX,NEW)
      RETURN
C
160   continue
      YPLOT=Y0PLOT-EARTHR*(PI-TH)*(DCOS(PH-PHR)*COSBET1
     1 +DSIN(PH-PHR)*SINBET1)                                      ! 15 Nov 2021
      XPLOT=X0PLOT-EARTHR*(PI-TH)*(DSIN(PH-PHR)*COSBET1
     1 -DCOS(PH-PHR)*SINBET1)                                      ! 15 Nov 2021
      CALL PLOT(XPLOT,YPLOT,NEW)                                   ! 11 Nov 2021
      RETURN                                                       ! 11 Nov 2021
C
10    SL=DSQRT(DMAX1(0.D0,1.D0-CL**2))
      TMP1=STH*dSIN (PH-PHL)
      TMP2=(CTH-CTHL*CL)/STHL
      ALPH2=0.D0
      IF (TMP1.NE.0.d0.OR.TMP2.NE.0.d0) ALPH2=DATAN2(TMP1,TMP2)
      CALL PLOT (EARTHR*CEA,EARTHR*dASIN(SL*dSIN (ALPH1-ALPH2)),NEW)
      RETURN
C
C          DRAW AXES AND CALL FOR LABELING AND TERMINATION OF THIS PLOT
C ----------------------------------------------------------------------!-------
      ENTRY ENDPLT
C ----------------------------------------------------------------------!-------
C
      GAMMA0 = R0/EARTHR                                          ! 29 Jun 2009
      IF(W(100) .EQ. 9.0d0) CALL CIRCLE(GAMMA0,LAM0,PH0)          ! 29 Jun 2009
      IF(W(100) .EQ. 9.0d0) CALL PLOTNOTE(W(101),IDSU)            ! 29 Jun 2009

      IF(W(125) .EQ. 3.0d0) CALL CIRCLE(GAMMA1,LAM1,PH1)          ! 29 Jun 2009
      IF(W(125) .EQ. 3.0d0) CALL PLOTNOTE(W(126),IDSDU)           ! 29 Jun 2009

      IF(W(100) .EQ. -6.0d0) CALL PLTNOTE(W(101),IDSU)            ! 14 Dec 2016
      IF(W(125) .EQ. -6.0d0) CALL PLTNOTE(W(126),IDSDU)           ! 14 Dec 2016

C     IF NEWWR IS STILL TRUE, NO PLOTS WHERE PRODUCED
      IF(NEWWR) RETURN
      NEWWR=.TRUE.
C
C     SIGNAL END OF PLOT
      IF(NPLTDP.GT.0) WRITE(NPLTDP) 4,(NWDSRK,I=2,NWDSRK)
C
      IF(PLT.EQ.0.d0) RETURN
C
      TICKX=1.d-2*(YT-YB)
      DTIC=TIC*EARTHR
      CALL SETXY(APLT,-ALPHA,RMIN,ALPHA,RMAX)
C
      IF(APLT.EQ.4.d0) GO TO 200
      IF(APLT.EQ.2.d0) GO TO 25
      IF(APLT.EQ.5.d0) GO TO 25                                    ! 12 Nov 2021
C
C     CURVLINEAR PROJECTIONS(=1,3)
      CALL ARCTIC(-ALPHA,ALPHA,RMIN,TICKX,PLTHLB)
      CALL ARCTIC(-ALPHA,ALPHA,RMAX,-TICKX,PLTHLB)
      GO TO 300
C
C     CARTESIAN PROJECTION PUT IN BOTTOM BOUNDARY
200   CALL TIKLINE(-ALPHA,RMIN,ALPHA,RMIN,TIC,-TICKX,PLTHLB)
C     PUT IN TOP BOUNDARY
      CALL TIKLINE(-ALPHA,RMAX,ALPHA,RMAX,TIC,TICKX,PLTHLB)
C
300   TIKL=2.d-2*ALPHA
      LTIC=TICV
CCCCC      IF(APLT.EQ.3.d0) LTIC=TICV*PFACTR  original
C
C.....new correction, May 1991
C
      if (aplt.eq.3.) ltic = ticv*pfactr
C     &               *(rmin/(earthr+(rmin-earthr)*pfactr))    ! 6 JULY 2004

C     PUT IN LEFT BOUNDARY
      CALL TIKLINE(-ALPHA,RMIN,-ALPHA,RMAX,LTIC,TIKL,PLTLB)
C     PUT IN RIGHT BOUNDARY
      CALL TIKLINE(ALPHA,RMIN,ALPHA,RMAX,LTIC,-TIKL,PLTLBN)
C
      GO TO 50
C   DRAW TICKS, BOX FOR HORIZONTAL PLOT
25    CALL DRAWTKS(DTIC,TICV,XL,XR,YB,YT,PLOT)
C
!50    IF(APLT.EQ.2.d0) THEN
50    IF(APLT.EQ.2.d0.or.APLT.EQ.5.d0) THEN                        ! 12 Nov 2021
        CALL PLTANOT(ID,OW/PIT2,XFQMDL,YFQMDL,XL,YB,XR,YT,DEGS,
     2     LLAT,LLON,RLAT,RLON,HB,HT,APLT,DTIC,DTIC/PFACTR,PLOT)
      ELSE
        CALL PLTANOT(ID,OW/PIT2,XFQMDL,YFQMDL,-ALPHA,RMIN,ALPHA,RMAX,
     2     DEGS,LLAT,LLON,RLAT,RLON,HB,HT,APLT,DTIC,TICV,PLOT)
      ENDIF
C
      CALL LABPLT
      CALL PLTEND
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE RAYPLT
C ----------------------------------------------------------------------!-------
 



C ----------------------------------------------------------------------!-------
      SUBROUTINE CIRCLE(gamma0,lam0,ph0)                         ! 29 Jun 2009
C ----------------------------------------------------------------------!-------
C Subroutine to draw a circle around the core of a vortex in a projection 
C on the ground
      implicit double precision (a-h,o-z)

      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10

C     COMMON DECK "FLAG" INSERTED HERE
      LOGICAL NEWWR,NEWWP,NEWTRC,PENET
      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH
      COMMON/FLGP/NSET

C     COMMON DECK "RKAM" INSERTED HERE
      DOUBLE PRECISION KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

      COMMON/WW/DUM(10),MAXW,W(1000)
      DOUBLE PRECISION MAXSTP,MAXERR,INTYP,LLAT,LLON
      EQUIVALENCE (EARTHR,W(1)),(RAY,W(2)),(XMTRH,W(3)),(TLAT,W(4)),
     1 (TLON,W(5)),(OW,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),
     2 (AZ1,W(10)),(AZBEG,W(11)),(AZEND,W(12)),(AZSTEP,W(13)),
     3 (BETA,W(14)),(ELBEG,W(15)),(ELEND,W(16)),(ELSTEP,W(17)),
     8   (RUNSUP,W(18)),(RCVRH,W(20)),
     4 (ONLY,W(21)),(HOP,W(22)),(MAXSTP,W(23)),(PLAT,W(24)),(PLON,W(25))
     5,(HMAX,W(26)),(RAYFNC,W(29)),(EXTINC,W(33)),
     6 (HMIN,W(27)),(RGMAX,W(28)),
     8 (INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),
     6 (STEP1,W(44)),(STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),
     7 (SKIP,W(71)),(RAYSET,W(72)),(PRTSRP,W(74)),(HITLET,W(75))
     9 ,(BINRAY,W(76)),(PAGLN,W(77)),(PLT,W(81)),(PFACTR,W(82)),
     1 (LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))
     2,(TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))

       DOUBLE PRECISION LAM0

      if(plt .ne. 2.0d0) return                                   ! 29 Jun 2009
      R=EARTHR + 100.0
      TH0 = PID2 - LAM0
      COSTH0=DCOS(TH0)
      SINTH0=DSIN(TH0)
      SINGAM=DSIN(GAMMA0)
      COSGAM=DCOS(GAMMA0)
      TH=TH0-GAMMA0
      PH=PH0
      NEWTRC=.TRUE.
      CALL RAYPLT
      DO 100 I=3, 360, 3
      ALPHA = I/DEGS
      COSTH = COSTH0*COSGAM+SINTH0*SINGAM*DCOS(ALPHA)
      TH=DACOS(COSTH)
      SINPHPH0 = SINGAM*DSIN(ALPHA)/DSIN(TH)
      PH = PH0 +DASIN(SINPHPH0)
      CALL RAYPLT
 100  CONTINUE
      RETURN
C ----------------------------------------------------------------------!-------
       END SUBROUTINE CIRCLE
C ----------------------------------------------------------------------!-------




C ----------------------------------------------------------------------!-------
      SUBROUTINE PLOTNOTE(WLABEL,IDSU)
C ----------------------------------------------------------------------!-------
C Subroutine to add notes to the plot in a projection on the ground

      implicit double precision (a-h,o-z)

      COMMON /DD/ INT,IOR,IT,IS,IC,ICC,IX,IY

      INTEGER IDSU(20)

      DOUBLE PRECISION WLABEL(24)

      CHARACTER NUMSTG*80

      CHARACTER*80 LABEL,CHID,id
      common/ident/id

      COMMON/WW/DUM(10),MAXW,W(1000)
      DOUBLE PRECISION MAXSTP,MAXERR,INTYP,LLAT,LLON
      EQUIVALENCE (EARTHR,W(1)),(RAY,W(2)),(XMTRH,W(3)),(TLAT,W(4)),
     1 (TLON,W(5)),(OW,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),
     2 (AZ1,W(10)),(AZBEG,W(11)),(AZEND,W(12)),(AZSTEP,W(13)),
     3 (BETA,W(14)),(ELBEG,W(15)),(ELEND,W(16)),(ELSTEP,W(17)),
     8   (RUNSUP,W(18)),(RCVRH,W(20)),
     4 (ONLY,W(21)),(HOP,W(22)),(MAXSTP,W(23)),(PLAT,W(24)),(PLON,W(25))
     5,(HMAX,W(26)),(RAYFNC,W(29)),(EXTINC,W(33)),
     6 (HMIN,W(27)),(RGMAX,W(28)),
     8 (INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),
     6 (STEP1,W(44)),(STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),
     7 (SKIP,W(71)),(RAYSET,W(72)),(PRTSRP,W(74)),(HITLET,W(75))
     9 ,(BINRAY,W(76)),(PAGLN,W(77)),(PLT,W(81)),(PFACTR,W(82)),
     1 (LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))
     2,(TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))

      if(plt .ne. 2.0d0) return                                   ! 29 Jun 2009

C Optional additional labeling of raypath plots
      IF(WLABEL(9) .EQ. 0.D0) RETURN                              ! 18 June 2009

C      WRITE(LABEL,'(a)') ID

      label = NUMSTG(idsu,20,'(20a4)')

      CHID=LABEL
      LABEL=CHID(1:20)
      IOR=0
      INT=1
      IT=0
      IS=3
      IX=WLABEL(10)
      IY=WLABEL(11)
      IF(IY.NE.0) CALL DDTEXT (2,LABEL)
      LABEL=CHID(21:40)
      IX=WLABEL(12)
      IY=WLABEL(13)
      IF(IY.NE.0) CALL DDTEXT (2,LABEL)
      IX=WLABEL(14)
      IY=WLABEL(15)
      IF(IY.NE.0) CALL DDTEXT (2,'De-focusing')
      IX=WLABEL(16)
      IY=WLABEL(17)
      IF(IY.NE.0) CALL DDTEXT (2,'----Focusing')
      IX=WLABEL(18)
      IY=WLABEL(19)
      IF(IY.NE.0) CALL DDTEXT (2,'----Virtual Source')
      IS=1                                                        ! 23 Jun 2009
      INT=0
      RETURN
C ----------------------------------------------------------------------!-------
       END SUBROUTINE PLOTNOTE
C ----------------------------------------------------------------------!-------



C ----------------------------------------------------------------------!-------
      SUBROUTINE RAYPLUS
C ----------------------------------------------------------------------!-------
!     Puts a plus sign on a rayplot
C ----------------------------------------------------------------------!-------
      COMMON /DD/ INT,IOR,IT,IS,IC,ICC,IX,IY                       ! 12 Apr 2011

         ixsave = ix
         iysave = iy
         ix = ix - 29
         iy = iy -6
         call ddtext(1,'+')
         ix = ixsave
         iy = iysave
      return
C ----------------------------------------------------------------------!-------
      END SUBROUTINE RAYPLUS
C ----------------------------------------------------------------------!-------


