C file PCPROF3.FOR
C ----------------------------------------------------------------------!-------
Contains PRFPLT, PRFSUP
C ----------------------------------------------------------------------!-------
!  9 Jul 2018 - Added diagnostic printing for w(73)=5.d0
!  9 Jul 2018 - Correct calculation for extraneous roots
!  6 Jul 2018 - Correct calculation of wavelength and vertical wavenumber
!  5 Jul 2018 - Added common block /RINK3/ to PRFSUP
!  2 Jul 2018 - Added print statements in calculation of wavelength
!  2 Jul 2018 - Added print statements in calculation of vertical wavenumber
! 29 Jun 2018 - Added print statements in calculation of wavelength
! 29 Jun 2018 - Added print statements in calculation of vertical wavenumber
! 28 Jun 2018 - Modified calculation of vertical wavenumber
! 27 Jun 2018 - Use propagate(i) in common block /RINK4/ to know if evanescent
! 27 Jun 2018 - Added common block /RINK4/ to PRFSUP
! 27 Jun 2018 - Calculate wavelength and vertical wavenumber only when valid
! 26 Jun 2018 - Calculate roots correctly
! 25 Jun 2018 - Calculate wavelength from SOL in common block /RINK/ 
! 25 Jun 2018 - Added common block /RINK/ to PRFSUP
! 22 May 2018 - Adjusted format on print
! 25 Aug 2011 - Changes recommended by Larry Baker
! 11 Apr 2011 - Modified calculation of vertical wavenumber in evanescent region
!  8 Apr 2011 - Added common block /rink2/
!  8 Apr 2011 - Corrected initialization for evanescent cases
!  8 Apr 2011 - Corrected formulas for wavelength and phase speed
!  7 Apr 2011 - Changed call FINIT to a function call
C 23 Jun 2009 - Changed the default character size for DDTEXT calls to 1.
C 19 Mar 2008 - Corrected lon to llon
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 24 May 2007 - Replaced ! with blank at begining of string in CALL ANNFIL
C 26 Mar 2007 - Added ! to begining of string in CALL ANNFIL
C  2 May 2006 - Removed ! from beginning of string in CALL ANNFIL
C 10 Apr 2006 - altered axis label for phase integral.
C  7 Apr 2006 - axis label for phase integral less that 40 characters.
C  6 Apr 2006 - altered axis label for phase integral
C  5 Apr 2006 - altered initialization for phase integral
C 31 Mar 2006 - Correction
C 31 Mar 2006 - Added validity of ray theory
C 21 Mar 2006 - Adjusting labels on graph
C 21 Mar 2006 - Correction
C 20 Mar 2006 - Added some common blocks for DDTEXT
C 17 Mar 2006 - Added printing of frequency, azimuth, beta, receiver ht.
C 17 Mar 2006 - Changed profile from transmitter height to receiver height
C 17 Mar 2006 - Changed profile from transmitter to left lon and lat
C 17 Mar 2006 - renamed entry points to phspeed and iphspeed
C 16 Mar 2006 - Added kzb to /RIN/
C 16 Mar 2006 - Changed kz2 to kz2a
C 14 Mar 2006 - Added phase speed to FUNCTION PRFSUP
C 14 Mar 2006 - Added wavelength to FUNCTION PRFSUP
C 14 Mar 2006 - Added vertical wavenumber to FUNCTION PRFSUP
C 14 Mar 2006 - Added kz2 to /RIN/
C  3 Mar 2006 - Changed units of Brunt-Vaisala frequency to Hz
C  3 Mar 2006 - Changed units of acoustic cutoff frequency to Hz
C  3 Mar 2006 - Correction
C  1 Mar 2006 - Added square of Brunt-Vaisala frequency to FUNCTION PRFSUP
C 28 Feb 2006 - Added call disper
C 27 Feb 2006 - Correction
C 24 Feb 2006 - Added argument to SPEEDV
C 24 Feb 2006 - Added acoustic cutoff frequency to FUNCTION PRFSUP
C 24 Feb 2006 - Added Brunt-Vaisala frequency to FUNCTION PRFSUP
C 31 August 2004 - Modified for M/S
C 4 August 2004 - Modified for M/S
C 25 June 2004 - Added 360.0 to WNANGV if necessary to make it positive.
C 25 June 2004 - Corrected output value for ENTRY IANGV
C 25 June 2004 - Shortened wind speed/component/direction character strings
C 25 June 2004 - Added '!' to end of wind direction character string
C 24 June 2004 - added wind/current direction to FUNCTION PRFSUP
C 2 Jun 2004 - Replaced (A8 by (1X,A8 in seven places
C 2,3 Jun 2004 - Added ! to beginning of string in CALL ANNFIL
C 2 Jun 2004 - Changed Km to km
C 2 Jun 2004 - Changed (DEGS K) to {K)
C 2 Jun 2004 - Added a RETURN statement to the end of SUBROUTINE PRFPLT
C ----------------------------------------------------------------------!-------
      SUBROUTINE PRFPLT(PLTMODE,FINIT,FUNX,FUNY,TLTS,NPTS,MAP,WORK
     1    ,FGEN,XPLT,YPLT,MARK)
C ----------------------------------------------------------------------!-------

C     'PRFPLT' IS A MULTI-FUNCTION GRAPHICAL SUPPORT ROUTINE TO AID IN
C     THE PRODUCTION OF ANNOTATED MULTI-TRACE X-Y PLOTS.
C     TO MORE EASILY CONTROL OPTIONS THE ARRAY 'PLTMODE' CONTAINS
C     A SET OF LOGICAL FLAGS WHICH WILL BE 'EXECUTED' IN THIS
C     PASS THROUGH 'PRFPLT'. THE FOLLOWING LIST OF FLAGS ARE
C     CURRENTLY IMPLEMENTED.
C
C      NFILL
C           CAUSES WORK ARRAY TO BE FILLED WITH VALUES OBTAINED BY CALL
C           TO EXTERNAL FUNCTIONS:
C             'FINIT' INITIALIZE FUNCTION.  FINIT HAS THE VALUE OF
C                     AN ENTRY POINT IN 'PRFSUP',E.G. ISPEDV.
C             'FUNX'  COMPUTE X-VALUES as for FINIT.
C             'FUNY'  COMPUTE Y-VALUES as for FINIT.
C      INISCX
C           CAUSES X SCALING ARRAY 'XPLT' TO BE INITIALIZED TO +-INFINI
C      INISCY
C           CAUSES Y SCALING ARRAY 'YPLT' TO BE INITIALIZED TO +-INFINI
C      NSCALE
C           CAUSES EXTREME X-Y VALUES TO BE DETERMINED SO THAT NO POINT
C           WILL BE OUT OF RANGE. MANUALLY SPECIFIED VALUES WILL NOT BE
C           OVERWRITTEN PROVIDED THEY ARE SUFFICIENT TO
C           INCLUDE ALL DATA VALUES.
C      INSFRAM
C           INSERT A FRAME IN THE GRAPHICS OUTPUT
C      NORMPLT
C           NORMAL PLOT. THIS WILL INCLUDE AUTOSCALING IF THE INCREMENT
C           VALUE USED FOR AXIS ANNOTATION IS GIVEN AS ZERO. IF EXTREME
C           VALUES ARE GIVEN BUT INCREMENT IS ZERO, AUTOSCALING WILL BE
C           USED TO DETERMINE THE INCREMENT AND SUBTICK INTERVAL
C           BUT WILL NOT DETERMINE EXTREME VALUES.
C      ENDFRAM
C           CAUSES A FINAL END FRAME TO BE OUTPUT TO GRAPHICS FILE.
C      ENDPLOT
C           COMPLETES OUTPUT TO GRAPHICS FILE, TERMINATING PLOTS.
C      NXSNEG, NYSNEG
C           FOR SUBROUTINE SUPNEG.  DETERMINES HOW MANY AXIS VALUES
C           ARE NEGATIVE.  ONLY FOR DISSPLA PACKAGE.
C
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY                    ! 20 Mar 2006
C      print *, 'IN=',IN,'IOR=',IOR,'IN=',IT,'IS=',IS,'IC=',IC ICC,IX,IY

C     COMMON DECK "CONST" INSERTED HERE                      ! 20 Mar 2006
      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 "WWR" INSERTED HERE                        ! 20 Mar 2006
      PARAMETER (NWARSZ=1000) 
      COMMON/WW/ dum(10), MAXW,W(NWARSZ)
      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))
C
      external finit                                               !  7 Apr 2011
      integer finit                                                !  7 Apr 2011
C ----------------------------------------------------------------------!-------
      character*80 id
      character*80 idd ! 17 Mar 2006
      common /ident/ id
      CHARACTER*(*) TLTS(4)
      INTEGER MAP(2)
      DOUBLE PRECISION WORK(NPTS,2),FGEN(5),XPLT(3),YPLT(3)
C
      LOGICAL PLTMODE(12), NOTSET
C
      LOGICAL PSTART,FSTART
C
      INTEGER NFILL, INISCX, INISCY, NSCALE, INPLT, INSFRAM, NORMPLT,
     &        ENDFRAM, ENDPLOT, NFULL, NXSNEG, NYSNEG

      PARAMETER (NFILL=1,INISCX=2,INISCY=3,NSCALE=4,INPLT=5,INSFRAM=6,
     &           NORMPLT=7,ENDFRAM=8,ENDPLOT=9,NFULL=10,NXSNEG=11,
     &           NYSNEG=12)

      DATA PSTART,FSTART/2*.TRUE./
      DATA BIG/1.E38/
C
C
C*****CHECK TO SEE IF ANY OF THE ELEMENTS OF "PLTMODE" ARE SET TO TRUE.
C
      NOTSET = .TRUE.
      DO 10 I=1,12
         IF (PLTMODE(I)) NOTSET = .FALSE.
10    CONTINUE
      IF (NOTSET) PLTMODE(NFULL) = .TRUE.
C
C.....Calculate the values and fill WORK array with the values to be plotted.
C     NOTE:  This is where the models are called i.e. FUNX and FUNY which
C            calculate the speeds.
C
      IF (PLTMODE(NFILL)) THEN
	 IF(NPTS.LT.2) STOP 'ERROR: Less than 2 data points.'
!  7 Apr 2011         CALL FINIT(FGEN(5),TLTS(1),TLTS(2),TLTS(4))
         itest = FINIT(FGEN(5),TLTS(1),TLTS(2),TLTS(4))            !  7 Apr 2011
         DD=(FGEN(2)-FGEN(1))/(NPTS-1)
         DO 100 I=0,NPTS-1
            F=FGEN(1)+DD*I
            WORK(I+1,1)=dINT(FUNX(F)*FGEN(3)*1.E5+.5)*1.E-5
            WORK(I+1,2)=FUNY(F)*FGEN(4) 
100         continue
      ENDIF
C
      IF (PLTMODE(INISCX)) THEN
         XPLT(1)=BIG
         XPLT(2)=-BIG
         XPLT(3)=0.0
      ENDIF
C
      IF (PLTMODE(INISCY)) THEN
         YPLT(1)=BIG
         YPLT(2)=-BIG
         YPLT(3)=0.0
      ENDIF
C
      XPN=XPLT(1)
      XPX=XPLT(2)
      YPN=YPLT(1)
      YPX=YPLT(2)
      IF (PLTMODE(NSCALE) .and. PLTMODE(iniscx)) THEN
         CALL VMNMX(WORK(1,MAP(1)),NPTS,XPLT(1),XPLT(2))
      endif
      if (PLTMODE(nscale) .and. PLTMODE(iniscy)) then
         CALL VMNMX(WORK(1,MAP(2)),NPTS,YPLT(1),YPLT(2))
      endif
C
C
ccccc	  IF (.NOT. PLTMODE(INSFRAM)) THEN
ccccc	     IF(XPN.GT.XPLT(1) .OR. XPX.LT.XPLT(2))
ccccc	1		 STOP 'ERROR:  Value in X array out of range'
ccccc	     IF(YPN.GT.YPLT(1) .OR. YPX.LT.YPLT(2))
ccccc	1		 STOP 'ERROR:  Value in Y array out of range'
ccccc	  ENDIF
C

      IF(PSTART) THEN
         CALL PLTINI
         PSTART=.FALSE.
      ENDIF
C
      IF (PLTMODE(INSFRAM)) THEN
         IF(.NOT.FSTART) CALL ENDPL(0)
C
C*****THE 100 IS NOT THE LENGTH OF THE CHAR STRING.  IT INDICATES THAT
C*****THE STRING IS TERMINATED BY A '$'.  JAH.
C
	 CALL TITLXY(TLTS(1),100,TLTS(2),100,TLTS(3),100,XAXIS,YAXIS)
         FSTART=.FALSE.
      ENDIF
C
      IF (PLTMODE(NORMPLT)) THEN
         IF(XPLT(3).LE.0.0) THEN
            CALL PSCALE(XPLT(1),XPLT(2),XPLT(3),NXTKS,XAXIS)

            IF(XPX.GT.XPN) THEN
                XPLT(1)=XPN
                XPLT(2)=XPX
            ENDIF
            CALL XTICKS(NXTKS)
         ENDIF
         IF(YPLT(3).LE.0.0) THEN
            CALL PSCALE(YPLT(1),YPLT(2),YPLT(3),NYTKS,YAXIS)
            IF(YPX.GT.YPN) THEN
                YPLT(1)=YPN
                YPLT(2)=YPX
            ENDIF
            CALL YTICKS(NYTKS)
         ENDIF
C
C.....Output the subheading before writing the data for the curve.
C     Subheading is ID in COMMON IDENT.  See READW for more.  It
C     contains additional user-supplied information about the particular plot.

      IOR=0           ! 20 Mar 2006
      IS=1                                                        ! 23 Jun 2009
      IX=0            ! 20 Mar 2006
      IY=0            ! 20 Mar 2006
      ix=400          ! 21 Mar 2006
      iy=1200         ! 20 Mar 2006
         call ddtext(8,id)

      freq=ow/pit2                                               ! 17 Mar 2006
      azdeg=az1*degs                                             ! 21 Mar 2006
      betadeg=beta*degs                                          ! 21 Mar 2006
      rechite=dabs(rcvrh)                                        ! 17 Mar 2006

      WRITE(idd,'(1X,F13.8," Hz,",f7.2," deg az,",f9.4," deg elev at "
     & ,f5.1," km")') freq,azdeg,betadeg,rechite   ! 17 Mar 2006 ! 22 May 2018

c      PRINT *,'IY = ',IY
      IY = IY + 300                                              ! 21 Mar 2006

c      PRINT *,'IY = ',IY

         call ddtext(8,idd)

         IF (PLTMODE(INSFRAM)) THEN
            NMBXS=1.5+(XPLT(2)-XPLT(1))/XPLT(3)
            NMBYS=1.5+(YPLT(2)-YPLT(1))/YPLT(3)
            NMBS=0
            IF (PLTMODE(NYSNEG)) THEN
               IDEL=NMBXS
               NMBS=NMBYS
            ENDIF
            IF (PLTMODE(NXSNEG)) THEN
               IDEL=0
               NMBS=NMBS+NMBXS
            ENDIF
            IF(NMBS.GT.0) CALL SUPNEG(IDEL,NMBS)
            CALL GRAFB(XPLT(1),XPLT(3),XPLT(2),YPLT(1),YPLT(3),YPLT(2))
         ENDIF
         MRK=0
         IF(MARK.LT.0) MRK=-1
C     DISSPLA CONVENTION IS TO AUTO-STEP THE MARKER SYMBOL.
C
         CALL CURVE(WORK(1,MAP(1)),WORK(1,MAP(2)),NPTS,MRK)
      ENDIF
C
      IF (PLTMODE(ENDFRAM) .AND. .NOT.FSTART) THEN
         CALL ENDPL(0)
         FSTART=.TRUE.
      ENDIF
C
      IF (PLTMODE(ENDPLOT) .AND. .NOT.PSTART) THEN
         CALL DONEPL
         PSTART=.TRUE.
      ENDIF
      RETURN ! 2 Jun 2004
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PRFPLT
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      INTEGER FUNCTION PRFSUP()
C ----------------------------------------------------------------------!-------
CONTAINS ENTRY POINTS TO CALCULATE THE PROFILE VALUES FOR EACH KIND OF PROFILE 
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

      TYPE :: value_and_derivs                                    ! 24 Feb 2006
         REAL(KIND=8) :: value
         REAL(KIND=8) :: pt,pr,pth,pph
      END TYPE 

      TYPE (value_and_derivs) :: kA2                       ! 24 Feb 2006
      TYPE (value_and_derivs) :: Ns                        ! 24 Feb 2006

      REAL(KIND=8) :: LPOLAR,LPOLRI,KPHK,KPHKI,KAY2,KAY2I 
     & ,kz2ar,kz2ai,kzbr,kzbi,kzr,kz2r                     ! 31 Mar 2006
c     & ,kz2a,kz2ai,kzb,kzbi,kz,kz2                        ! 31 Mar 2006
      COMPLEX*16 kz2a, kzb, kz, kz2                        ! 31 Mar 2006

      COMMON/RINPL/dispm                                   ! 24 Feb 2006
      CHARACTER(LEN=6) :: dispm                            ! 24 Feb 2006
      COMMON/RIN/modrin,rayname(2,3),type(3),space         ! 24 Feb 2006
      COMMON/RIN/OMEGMIN,OMEGMAX,KAY2,KAY2I 
     & ,kz2ar,kz2ai,kzbr,kzbi                              ! 31 Mar 2006
     1 ,H,HI,PHT,PHTI,PHR,PHRI,PHTH,PHTHI,PHPH,PHPHI       ! 14 Mar 2006
     2 ,PHOW,PHOWI,PHKR,PHKRI,PHKTH,PHKTI,PHKPH,PHKPI      ! 14 Mar 2006
     3 ,KPHK,KPHKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN          ! 14 Mar 2006
c     & , kz2a,kzb                                         ! 31 Mar 2006
      EQUIVALENCE (kz2a,kz2ar),(kzb,kzbr),(kz,kzr), (kz2,kz2r)     ! 31 Mar 2006

      COMMON /rink2/ kr2,kth2,kph2                                 !  8 Apr 2011

      CHARACTER(LEN=64) :: modrin                          ! 24 Feb 2006
      LOGICAL :: space                                     ! 24 Feb 2006

! 25 Jun 2018 - Added common block /RINK/ to PRFSUP

      common /rink/ kaysq(4),kaysqb(4),kaysqa(4)                   ! 14 Jun 2010
     & ,kaysqi(4),kaysqbi(4),kaysqai(4)                            !  5 Aug 2010
     & ,sol(4),soli(4),solb(4),solbi(4),sola(4),solai(4)           ! 15 Jun 2010
     & ,alpha,dd                                                   ! 28 Jul 2010
     & ,relerr(4),relerra(4),relerrb(4)                            ! 28 Jun 2010
     & ,errRell(4),errRela(4),errRelb(4)                           ! 19 Nov 2010
     & ,n3,n4,nbig,nsmall                                          ! 15 Jul 2010
     & ,Hrel,errRel                                                ! 27 Jul 2010
      complex*16 relerr,relerra,relerrb                            ! 28 Jun 2010
      complex*16 errRell,errRela,errRelb                           ! 19 Nov 2010
      double precision kaysq,kaysqb,kaysqa                         ! 14 Jun 2010
     & ,kaysqi,kaysqbi,kaysqai                                     !  5 Aug 2010
      dimension dd(0:4)
     & ,ddz(0:4)                                                   ! 11 Jul 2017
!     & ,sol(4),soli(4),solb(4),solbi(4),sola(4),solai(4)          ! 15 Jun 2010

! 27 Jun 2018 - Added common block /RINK4/ to PRFSUP

      integer propagate(4)                                         ! 27 Jun 2018
      COMMON /rink4/ propagate                                     ! 27 Jun 2018

!  5 Jul 2018 - Added common block /RINK3/ to PRFSUP

      DOUBLE PRECISION azz,bzz,czz                                 !  5 Aug 2016
      COMMON /rink3/ azz,bzz,czz                                   !  5 Aug 2016
     &,czz0a,czz0b,czzf0,czzfth,czzfph,czzfthth,czzfphph,czzfthph  !  5 Aug 2016

C     COMMON DECK "UU" INSERTED HERE
      DOUBLE PRECISION MODU
      COMMON/UU/MODU(4)
     1  ,V  ,PVT  ,PVR  ,PVTH  ,PVPH
     2  ,VR ,PVRT ,PVRR ,PVRTH ,PVRPH
     3  ,VTH,PVTHT,PVTHR,PVTHTH,PVTHPH
     4  ,VPH,PVPHT,PVPHR,PVPHTH,PVPHPH

C     COMMON DECK "B1" INSERTED HERE
      INTEGER UMX,UNTBL,UITBL,UFRMTBL,IDSU(10)
      COMMON/B1/UMX,UNTBL(10),UITBL(10),UFRMTBL(10),UGP(10)
      EQUIVALENCE (UGP,IDSU)

C     COMMON DECK "CC" INSERTED HERE
      DOUBLE PRECISION MODC
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH

C     COMMON DECK "B3" INSERTED HERE
      INTEGER CMX,CNTBL,CITBL,CFRMTBL,IDSC(10)
      COMMON/B3/CMX,CNTBL(10),CITBL(10),CFRMTBL(10),CGP(512)
      EQUIVALENCE (CGP,IDSC),(ANC,CGP(11))

C     COMMON DECK "TT" INSERTED HERE
      DOUBLE PRECISION MODT
      COMMON/TT/MODT(4), T,PTT,PTR,PTTH,PTPH

C     COMMON DECK "B5" INSERTED HERE
      INTEGER TMX,TNTBL,TITBL,TFRMTBL,IDST(10)
      COMMON/B5/TMX,TNTBL(10),TITBL(10),TFRMTBL(10),TGP(262)
      EQUIVALENCE (TGP,IDST),(ANT,TGP(11))

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 "GG" INSERTED HERE
      DOUBLE PRECISION MODG
      COMMON/GG/MODG(4)
      COMMON/GG/G,PGR,PGRR,PGRTH,PGRPH
      COMMON/GG/PGTH,PGPH,PGTHTH,PGPHPH,PGTHPH,GSELECT,GTIME

C     COMMON DECK "RKAM" INSERTED HERE

      DOUBLE PRECISION KR,KTH,KPH,kth2,kph2                        ! 14 Mar 2006
     &               , kr2                                         !  8 Apr 2011
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

C     COMMON DECK "WWR" INSERTED HERE
      PARAMETER (NWARSZ=1000) 
      COMMON/WW/ dum(10), MAXW,W(NWARSZ)
      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))
     3,(SUBINT,W(95))
C
C     COMMON DECK "CPROCFL" INSERTED HERE
      INTEGER PMX,PNTBL,PITBL,PFRMTBL,IDSP(10)
C     PARAMETER DECK "PGROUPS"
      PARAMETER (NCHPG1=11,NWPV=500,NSPGP=NCHPG1+2*NWPV+1)         ! 25 Aug 2011
      PARAMETER (MNGRP=9,MXGRP=69,MXLIST=MXGRP-MNGRP+2)
      COMMON/PROCFL/LIST(MXLIST)
      COMMON/PROCFL/PMX,PNTBL(10),PITBL(10),PFRMTBL(10),PGP(NSPGP)
      EQUIVALENCE (PGP,IDSP)
C
      CHARACTER ITOC*7
      CHARACTER*40 TLTP,TLTY,TLTL
      INTEGER STRIM
C
      CALL ITOPOG
      R=EARTHR
      TH=PID2-LLAT
      PH=LLON
      PRFSUP=0
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY ISPEDV(X,TLTP,TLTY,TLTL)
      CALL ISPEED
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODC
! 2 May 2006     CALL ANNFIL('!S#OUND !S#PEED ('//TLTY(:STRIM(TLTY))//'/SEC!)' 
! 2 May 2006     1         ,TLTY) 
      CALL ANNFIL('S#OUND !S#PEED ('//TLTY(:STRIM(TLTY))//'/SEC!)'  ! 2 May 2006
     1         ,TLTY)                                              !  2 May 2006
C*****CHANGED 10a8 TO 10a4.
      WRITE(TLTL,'(10A4)') (CGP(I),I=1,10)
      ISPEDV=1
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY SPEEDV(x)                             ! added argument - 24 Feb 2006
      r=x                                                  ! added - 24 Feb 2006
      CALL SPEED
      SPEEDV=SQRT(CS)
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY IWINDV(X,TLTP,TLTY,TLTL)
      FX=1.D0                        ! 31 August 2004
      IF(TLTY(1:1).EQ.'M')FX=1.D3    ! 31 August 2004
      CALL IWINDR
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODU
c      CALL ANNFIL('!W#IND !S#PEED (' ! 25 June 2004 ! 2 May 2006
c     1//TLTY(:STRIM(TLTY))//'/SEC!)',TLTY) ! 2 May 2006
      CALL ANNFIL('W#IND !S#PEED (' ! 25 June 2004 ! 2 May 2006
     1//TLTY(:STRIM(TLTY))//'/SEC!)',TLTY) ! 2 May 2006
C*****CHANGED 10A8 TO 10A4
      WRITE(TLTL,'(10A4)') (UGP(I),I=1,10)
      IWINDV=1
      RETURN

C ----------------------------------------------------------------------!-------
      ENTRY IEASTV(X,TLTP,TLTY,TLTL)
      FX=1.D0                        ! 4 August 2004
      IF(TLTY(1:1).EQ.'M')FX=1.D3    ! 4 August 2004
      CALL IWINDR
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODU
C*****CHANGED 10A8 TO 10A4
! 2 May 2006      call ANNFIL('!E#ASTWARD !W#IND !C#OMP. (' ! 25 June 2004 
! 2 May 2006     &		  //TLTY(:STRIM(TLTY))//'/SEC!)' 
! 2 May 2006     &		   ,TLTY) 
      call ANNFIL('E#ASTWARD !W#IND !C#OMP. ('       ! 25 June 2004 ! 2 May 2006
     &		  //TLTY(:STRIM(TLTY))//'/SEC!)'                  ! 2 May 2006
     &		   ,TLTY)                                         ! 2 May 2006
      WRITE(TLTL,'(10A4)') (UGP(I),I=1,10)
      IEASTV=1
      RETURN

C ----------------------------------------------------------------------!-------
      ENTRY ISOUTV(X,TLTP,TLTY,TLTL)
      FX=1.D0                        ! 31 August 2004
      IF(TLTY(1:1).EQ.'M')FX=1.D3    ! 31 August 2004
      CALL IWINDR
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODU
C*****CHANGED 10A8 TO 10A4
! 2 May 2006      call ANNFIL ('!S#OUTHWARD !W#IND !C#OMP. (' ! 25 June 2004 
! 2 May 2006     &		   //TLTY(:STRIM(TLTY))//'/SEC!)' 
! 2 May 2006     &		   ,TLTY) 
      call ANNFIL ('S#OUTHWARD !W#IND !C#OMP. ('     ! 25 June 2004 ! 2 May 2006
     &		   //TLTY(:STRIM(TLTY))//'/SEC!)'                 ! 2 May 2006
     &		   ,TLTY)                                         ! 2 May 2006
      WRITE(TLTL,'(10A4)') (UGP(I),I=1,10)
      ISOUTV = 1
      RETURN

C ----------------------------------------------------------------------!-------
      ENTRY IVERTV(X,TLTP,TLTY,TLTL)
      FX=1.D0                                                   ! 31 August 2004
      IF(TLTY(1:1).EQ.'M')FX=1.D3                               ! 31 August 2004
      CALL IWINDR
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODU
C*****CHANGED 10A8 TO 10A4
c      call ANNFIL('!V#ERTICAL !W#IND !C#OMP. (' ! 25 June 2004 ! 2 May 2006
c     &		  //TLTY(:STRIM(TLTY))//'/SEC!)' ! 2 May 2006
c     &		   ,TLTY) ! 2 May 2006
      call ANNFIL('V#ERTICAL !W#IND !C#OMP. (' ! 25 June 2004 ! 2 May 2006
     &		  //TLTY(:STRIM(TLTY))//'/SEC!)' ! 2 May 2006
     &		   ,TLTY) ! 2 May 2006
      WRITE(TLTL,'(10A4)') (UGP(I),I=1,10)
      IVERTV = 1
      RETURN

C ----------------------------------------------------------------------!-------
      ENTRY IANGV(X,TLTP,TLTY,TLTL) ! 24 June 2004
      CALL IWINDR
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODU
C*****CHANGED 10A8 TO 10A4
! 2 May 2006      call ANNFIL('!W#IND !D#IRECTION (!D#EG.!)'      ! 25 June 2004 
! 2 May 2006     &		   ,TLTY)                           ! 24 June 2004 
      call ANNFIL('W#IND !D#IRECTION (!D#EG.!)'      ! 25 June 2004 ! 2 May 2006
     &		   ,TLTY)                          ! 24 June 2004 ! 2 May 2006
      WRITE(TLTL,'(10A4)') (UGP(I),I=1,10)
      IANGV = 1                                                   ! 25 June 2004
      RETURN

C
C Wind (or current) magnitude.

C ----------------------------------------------------------------------!-------
      ENTRY WNMAGV
      CALL WINDR
      WNMAGV=SQRT(VR*VR+VTH*VTH+VPH*VPH)*FX                     ! 31 August 2004
      RETURN

C Wind (or current) vertical component.

C ----------------------------------------------------------------------!-------
      ENTRY WNVRTV
      CALL WINDR
      WNVRTV=VR*FX                                              ! 31 August 2004
      RETURN

C Wind (or current) southward component.

C ----------------------------------------------------------------------!-------
      ENTRY WNSTHV
      CALL WINDR
      WNSTHV=VTH*FX                                             ! 31 August 2004
      RETURN

C Wind (or current) eastward component.

C ----------------------------------------------------------------------!-------
      ENTRY WNESTV
      CALL WINDR
      WNESTV=VPH*FX                                              ! 4 August 2004
      RETURN

C Wind (or current) direction.                                      24 June 2004

C ----------------------------------------------------------------------!-------
      ENTRY WNANGV                                                ! 24 June 2004
      CALL WINDR
      WNANGV=DEGS*ATAN2(-VPH,VTH)                                 ! 24 June 2004
      IF(WNANGV.LT.0.D0) WNANGV = WNANGV + 360.0D0                ! 25 June 2004
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY COLAT
      COLAT=PID2-TH
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY REARTH(X)
      R=X
      CALL TOPOG
      REARTH=R-EARTHR
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY ITEMPV(X,TLTP,TLTY,TLTL)
      CALL ITEMP
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODT
! 2 May 2006      CALL ANNFIL('!T#EMPERATURE (!K)',TLTY)                    
! 2 May 2006      CALL ANNFIL('T#EMPERATURE (!K)',TLTY)   
! 26 Mar 2006      CALL ANNFIL('!T#EMPERATURE (!K)',TLTY)   
       CALL ANNFIL(' T#EMPERATURE (!K)',TLTY)                      ! 24 May 2007
C*****CHANGED 10A8 TO 10A4
      WRITE(TLTL,'(10A4)') (TGP(I),I=1,10)
      ITEMPV=1
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY TOPV(X)
      TH=X
      CALL TOPOG
      TOPV=-G/PGR
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY TEMPV(X)
      R=X
      CALL TEMP
      TEMPV=T
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY IPTPLT(X,TLTP,TLTY,TLTL)
      NP=X-(MNGRP-2)
C
C     INITIALIZE OFFSETS FOR X AND Y ARRAYS
C     NOTE:  If error message gives number of points = -7, that means that
C            X must be 0.0.

      IF(NP.LT.2 .OR. NP.GT.LIST(1)) then
     	   PRINT*,'ERROR: Invalid number of points in tabular data: ',np
           STOP
      endif
C
      TLTP='DATA POINT MODEL'//ITOC(NP+(MNGRP-2))
C
      NLN=LIST(NP)

      IF(NLN.LE.0) STOP 'ERROR: Invalid LIST (see PRFSUP routine).'
C
C*****NOTE:  THE LOOP FOR PGP IS 1 TO 10 SINCE NCHPG1=11 AND NLN=1
      WRITE(TLTL,'(10A4)') (PGP(I),I=NLN,NLN+(NCHPG1-2))
C
C     GET SET POINTER AND STEP OVER IDENTIFIER STRING
      NLN=NLN+(NCHPG1-1)
C
C     GET NUMBER OF AVAILABLE POINTS
      NPTS=PGP(NLN)
      N2=NPTS/2
      IPTPLT=N2
C     X VALUES ARE NEXT
      NLX=NLN+1
C     THEN Y VALUES
      NLY=NLX+N2
C     SET LAST AVAILABLE LIMIT
      NMX=NLY+N2
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY PTPLTX(X)
      N=NLX+X
      IF(N.LT.NLX .OR. N.GE.NLY)
     &	 STOP 'ERROR: Invalid index into X array'
      PTPLTX=PGP(N)
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY PTPLTY(X)
      N=NLY+X
      IF(N.LT.NLY .OR. N.GT.NMX)
     &	 STOP 'ERROR: Invalid index into Y array'
      PTPLTY=PGP(N)
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY ITOPV (X,TLTP,TLTY,TLTL)
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODG
      tlty = '                         '
! 2 May 2006      CALL ANNFIL('!R#ANGE !(#KM!)',TLTY) 
      CALL ANNFIL('R#ANGE !(#KM!)',TLTY)                            ! 2 May 2006
      ITOPV=1
      RETURN

C ----------------------------------------------------------------------!-------
C Brunt-Vaisala frequency -                                          24 Feb 2006
      ENTRY IBV(X,TLTP,TLTY,TLTL)
      CALL IDISPER
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODT
      WRITE(TLTL,'(10A4)') (TGP(I),I=1,10)
c      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') dispm
! 2 May 2006      CALL ANNFIL('!B#RUNT-!V#AISALA FREQUENCY ('// 
! 2 May 2006     1'!H#Z!)'         ,TLTY)                           ! 3 Mar 2006
      CALL ANNFIL('B#RUNT-!V#AISALA FREQUENCY ('//                  ! 2 May 2006
     1'!H#Z!)'         ,TLTY)                          ! 3 Mar 2006 ! 2 May 2006
c      WRITE(TLTL,'(A40)') modrin
      IBV=1
      RETURN

C ----------------------------------------------------------------------!-------
C Brunt-Vaisala frequency - 24 Feb 2006
      ENTRY BV(x)
      r=x
c      print *, ' r = ', r
      call disper                              ! 28 Feb 2006
      CALL Brunt(Ns)
c      print *, 'Ns = ', Ns%value
c      BV=Ns%value
      IF(Ns%value.GE.0.) BV=DSQRT(Ns%value)/PIT2                  ! 3 Mar 2006
      IF(Ns%value.LT.0.) BV=-DSQRT(-Ns%value)/PIT2                ! 3 Mar 2006
c      print *, 'bv = ', bv
      RETURN

C ----------------------------------------------------------------------!-------
C Square of Brunt-Vaisala frequency -  1 Mar 2006
      ENTRY IBV2(X,TLTP,TLTY,TLTL)
      CALL IDISPER
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODT
      WRITE(TLTL,'(10A4)') (TGP(I),I=1,10)
c      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') dispm
c      CALL ANNFIL('!B#RUNT-!V#AISALA FREQUENCY ('//
c     1'!H#Z!)'         ,TLTY)                                      ! 3 Mar 2006
! 2 May 2006      CALL ANNFIL('!N#SQUARED ('// 
! 2 May 2006     1'!H#Z !H#Z!)'         ,TLTY)                      ! 3 Mar 2006 
      CALL ANNFIL('N#SQUARED ('// ! 2 May 2006
     1'!H#Z !H#Z!)'         ,TLTY)                     ! 3 Mar 2006 ! 2 May 2006
c      WRITE(TLTL,'(A40)') modrin
      IBV2=1                                             ! correction 3 Mar 2006
      RETURN

C ----------------------------------------------------------------------!-------
C Square of Brunt-Vaisala frequency -                                 1 Mar 2006
      ENTRY BV2(x)
      r=x
c      print *, ' r = ', r
      call disper                                                  ! 28 Feb 2006
      CALL Brunt(Ns)
c      print *, 'Ns = ', Ns%value
      BV2=Ns%value/(PIT2*PIT2)                           ! correction 3 Mar 2006
c      print *, 'bv = ', bv
      RETURN

C ----------------------------------------------------------------------!-------
C Acoustic cutoff frequency -                                        24 Feb 2006
      ENTRY IACF(X,TLTP,TLTY,TLTL)
      CALL IDISPER
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODT
      WRITE(TLTL,'(10A4)') (TGP(I),I=1,10)
c      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') dispm
! 2 May 2006      CALL ANNFIL('!A#COUSTIC !C#UTOFF FREQUENCY ('// 
! 2 May 2006     1'!H#Z!)'         ,TLTY)                           ! 3 Mar 2006
      CALL ANNFIL('A#COUSTIC !C#UTOFF FREQUENCY ('//                ! 2 May 2006
     1'!H#Z!)'         ,TLTY)                          ! 3 Mar 2006 ! 2 May 2006
c      WRITE(TLTL,'(A40)') modrin
      IACF=1
      RETURN

C ----------------------------------------------------------------------!-------
C Acoustic cutoff frequency -                                        24 Feb 2006
      ENTRY ACF(x)
      r=x
c      print *, ' r = ', r
      call disper                                                  ! 28 Feb 2006
      CALL speed
c      print *, ' Cs = ', Cs
      CALL cutoff(kA2)
c      print *, 'kA2 profile = ', kA2%value
c      print *,'h profile = ',r
      ACF=DSQRT(Cs*kA2%value)/PIT2                                  ! 3 Mar 2006
c      print *, 'acf = ', acf
      RETURN
C
C ----------------------------------------------------------------------!-------
C Vertical wavenumber - 14 Mar 2006
      ENTRY IVWN(X,TLTP,TLTY,TLTL)
      CALL IDISPER
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODT
      WRITE(TLTL,'(10A4)') (TGP(I),I=1,10)
! 2 May 2006      CALL ANNFIL('!V#ERTICAL !W#AVENUMBER ('// 
! 2 May 2006     1'#/KM!)'         ,TLTY) 
      CALL ANNFIL('V#ERTICAL !W#AVENUMBER ('//                      ! 2 May 2006
     1'#/KM!)'         ,TLTY)                                       ! 2 May 2006
c      WRITE(TLTL,'(A40)') modrin

      r=earthr + dabs(rcvrh)                                       ! 17 Mar 2006
      th=pid2-llat                                                 ! 17 Mar 2006
      ph=llon                                                      ! 17 Mar 2006
      kr=dsin(beta)
      kth=-dsin(pid2-beta)*dcos(az1)
      kph= dsin(pid2-beta)*dsin(az1)
      call disper
      kr2=kr*kr*abs(kay2)                            !  8 Apr 2011 ! 27 Jun 2018
      kth2=kth*kth*abs(kay2)                         !  8 Apr 2011 ! 27 Jun 2018
      kph2=kph*kph*abs(kay2)                         !  8 Apr 2011 ! 27 Jun 2018
      IVWN=0
      IIVWN=IVWN                                                   ! 27 Jun 2018
      if(ray.eq.0.d0.and.propagate(n3).eq.0) return                ! 27 Jun 2018
      if(ray.eq.-1.d0.and.propagate(nsmall).eq.0) return           ! 27 Jun 2018
      if(ray.eq.-2.d0.and.propagate(nbig).eq.0) return             ! 27 Jun 2018
      if(ray.ge.1.d0.and.ray.le.4.d0.and.propagate(int(ray)).eq.0) 
     & return                                                      ! 27 Jun 2018
c      call renorm(kr,kay2,3)
      kr=dsqrt(abs(kay2))*dsin(beta)                               ! 27 Jun 2018
      up=sign(1.d0,kr)                                             !  6 Jul 2018
      kth=-dsqrt(abs(kay2))*dsin(pid2-beta)*dcos(az1)              ! 27 Jun 2018
      kph= dsqrt(abs(kay2))*dsin(pid2-beta)*dsin(az1)              ! 27 Jun 2018
      if(W(73).eq.5.d0) print*,'wavenumber root',int(ray),
     & 'KR = ',KR, 'KTH = ',KTH, 'KPH = ',KPH                      !  2 Jul 2018
      if(ray.ge.1.d0.and.ray.le.4.d0.and.W(73).eq.5.d0) print*,
     & 'sol= ',sol(int(ray)),soli(int(ray))                        !  9 Jul 2018
!     & ,'r=',r,'th=',th,'ph=',ph,'kay2=',kay2                      ! 17 Jul 2018
      if(kay2.le.0.d0) then                                        !  9 Jul 2018
      kr=-kr
      kth=-kth
      kph=-kph
      up=-up
      endif                                                        !  9 Jul 2018
      IVWN=1
      IIVWN=IVWN                                                   ! 27 Jun 2018
      RETURN

C ----------------------------------------------------------------------!-------
C Vertical wavenumber                                              - 14 Mar 2006
      ENTRY VWN(x)
      VWN=0.D0                                                     ! 27 Jun 2018
      IF(IIVWN.EQ.0) RETURN                                        ! 27 Jun 2018
      r=x
      if(RCVRH.lt.0.d0) print*,'RCVRH = ',RCVRH
      if(kzbr.ne.0.d0) print*,'kzbr before = ',kzbr

      call disper
      if(kzbr.ne.0.d0) print*,'kzbr after = ',kzbr
c      if(VWN.lt.0.d0) print*,'VWN second = ',VWN

      VWN = 10.D0                                                  ! 28 Jun 2018

      if(azz.eq.0.d0.and.bzz.eq.0.d0) return                       !  6 Jul 2018
      if(azz.eq.0.d0) then                                         !  6 Jul 2018
      kr = -czz/(2.d0*bzz)                                         !  6 Jul 2018
      vwn=kr
      return
      endif                                                        !  6 Jul 2018
      rad = bzz**2 + azz*czz                                       !  6 Jul 2018
      if(rad.ge.0.d0) then                                         !  6 Jul 2018
      kr = (bzz+up*sqrt(rad))/azz 
      else
      kr = 5.D-3                                                
      endif                                                        !  6 Jul 2018
      vwn=kr
      
      height=r-earthr                                              !  2 Jul 2018
      if(W(73).eq.5.d0) print*,'height=',height,'kay2=',kay2
     & ,'kz2ar=',kz2ar,'KR = ',KR, 'KTH = ',KTH, 'KPH = ',KPH      !  9 Jul 2018

      RETURN
C
C ----------------------------------------------------------------------!-------
C wavelength - 14 Mar 2006
      ENTRY IWL(X,TLTP,TLTY,TLTL)
c      PRINT*,'IWL REACHED'
      CALL IDISPER
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODT
      WRITE(TLTL,'(10A4)') (TGP(I),I=1,10)
! 2 May 2006      CALL ANNFIL('!W#AVELENGTH ('// 
! 2 May 2006     1'#KM!)'         ,TLTY) 
      CALL ANNFIL('W#AVELENGTH ('//                                 ! 2 May 2006
     1'#KM!)'         ,TLTY)                                        ! 2 May 2006
c      WRITE(TLTL,'(A40)') modrin

c      r=earthr+xmtrh                                              ! 31 Mar 2006
      r=earthr + dabs(rcvrh)
c      print*,'r profile = ',r
      th=pid2-llat                                                 ! 17 Mar 2006
      ph=llon                                        ! 17 Mar 2006 ! 19 Mar 2008
c      print*,'profile beta = ', beta*degs
      kr=dsin(beta)
      kth=-dsin(pid2-beta)*dcos(az1)
      kph= dsin(pid2-beta)*dsin(az1)
c      print*,'profile kr = ',kr
c      print*,'profile kth = ',kth
c      print*,'profile kph = ',kph
      call disper
      kr2=kr*kr*abs(kay2)                            !  8 Apr 2011 ! 27 Jun 2018
      kth2=kth*kth*abs(kay2)                         !  8 Apr 2011 ! 27 Jun 2018
      kph2=kph*kph*abs(kay2)                         !  8 Apr 2011 ! 27 Jun 2018
      IWL=0
      IIWL=IWL                                                     ! 27 Jun 2018
      if(ray.eq.0.d0.and.propagate(n3).eq.0) return                ! 27 Jun 2018
      if(ray.eq.-1.d0.and.propagate(nsmall).eq.0) return           ! 27 Jun 2018
      if(ray.eq.-2.d0.and.propagate(nbig).eq.0) return             ! 27 Jun 2018
      if(ray.ge.1.d0.and.ray.le.4.d0.and.propagate(int(ray)).eq.0) 
     & return                                                      ! 27 Jun 2018
c      call renorm(kr,kay2,3)
c      PRINT*,'kay2 = ',kay2
      kr=dsqrt(abs(kay2))*dsin(beta)                               ! 27 Jun 2018
      up=sign(1.d0,kr)                                             !  6 Jul 2018
      kth=-dsqrt(abs(kay2))*dsin(pid2-beta)*dcos(az1)              ! 27 Jun 2018
      kph= dsqrt(abs(kay2))*dsin(pid2-beta)*dsin(az1)              ! 27 Jun 2018
      if(W(73).eq.5.d0) print*,'wavelength root',int(ray),
     & 'KR = ',KR, 'KTH = ',KTH, 'KPH = ',KPH                      !  2 Jul 2018
      if(ray.ge.1.d0.and.ray.le.4.d0.and.W(73).eq.5.d0) print*,
     & 'sol= ',sol(int(ray)),soli(int(ray))                        !  9 Jul 2018
!     & ,'r=',r,'th=',th,'ph=',ph,'kay2=',kay2                      ! 17 Jul 2018
      if(kay2.le.0.d0) then                                        !  9 Jul 2018
      kr=-kr
      kth=-kth
      kph=-kph
      up=-up
      endif                                                        !  9 Jul 2018
      IWL=1
      IIWL=IWL                                                     ! 27 Jun 2018
      RETURN

C ----------------------------------------------------------------------!-------
C wavelength - 14 Mar 2006
      ENTRY WL(x)
      WL=0.D0                                                      ! 27 Jun 2018
      IF(IIWL.EQ.0) RETURN                                         ! 27 Jun 2018
      r=x
      call disper
c      print*,'wavelength kz2ar = ',kz2ar
c      print*,'wavelength kz2r = ',kz2r

      kr = 10.d0                                                   ! 26 Jun 2018
      kr2 = kr*kr            
      TEMPk2=kr2+kth2+kph2                                         !  8 Apr 2011
c      print*,'wavelength TEMPk2 = ', TEMPk2
      WL=PIT2/DSQRT(TEMPk2)                       !  8 Apr 2011
 
      if(azz.eq.0.d0.and.bzz.eq.0.d0) return                       !  6 Jul 2018
      if(azz.eq.0.d0) then                                         !  6 Jul 2018
      kr = -czz/(2.d0*bzz)                                         !  6 Jul 2018
      kr2 = kr*kr            
      TEMPk2=kr2+kth2+kph2                                         !  8 Apr 2011
      WL=PIT2/DSQRT(TEMPk2)
      return
      endif                                                        !  6 Jul 2018
      rad = bzz**2 + azz*czz                                       !  6 Jul 2018
      if(rad.ge.0.d0) then                                         !  6 Jul 2018
      kr = (bzz+up*sqrt(rad))/azz 
      else
      kr = .005D0                                                
      endif                                                        !  6 Jul 2018
      kr2 = kr*kr            

      TEMPk2=kr2+kth2+kph2                                        !  8 Apr 2011
c      print*,'wavelength TEMPk2 = ', TEMPk2
      WL=PIT2/DSQRT(TEMPk2)                       !  8 Apr 2011

      height=r-earthr                                              !  2 Jul 2018
      if(W(73).eq.5.d0) print*,'height=',height,'kay2=',kay2
     & ,'kz2ar=',kz2ar,'KR = ',KR, 'KTH = ',KTH, 'KPH = ',KPH      !  9 Jul 2018

      RETURN

C ----------------------------------------------------------------------!-------
C phase speed - 14 Mar 2006
      ENTRY IPHSPEED(X,TLTP,TLTY,TLTL) ! 17 Mar 2006
      CALL IDISPER
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODT
      WRITE(TLTL,'(10A4)') (TGP(I),I=1,10)
c      CALL ANNFIL('!P#HASE SPEED ('// ! 2 May 2006
c     1'#KM/S!)'         ,TLTY) ! 2 May 2006
      CALL ANNFIL('P#HASE SPEED ('// ! 2 May 2006
     1'#KM/S!)'         ,TLTY) ! 2 May 2006
c      WRITE(TLTL,'(A40)') modrin

c      r=earthr+xmtrh                                              ! 31 Mar 2006
      r=earthr + dabs(rcvrh)
      th=pid2-llat                                                 ! 17 Mar 2006
      ph=llon                                                      ! 17 Mar 2006
      kr=dsin(beta)
      kth=-dsin(pid2-beta)*dcos(az1)
      kph= dsin(pid2-beta)*dsin(az1)
      call disper
      kr2=kr*kr*kay2                                               !  8 Apr 2011
      kth2=kth*kth*kay2                                            !  8 Apr 2011
      kph2=kph*kph*kay2                                            !  8 Apr 2011
      IPHSPEED=0     ! 17 Mar 2006
      if(kay2.le.0.d0) return
c      call renorm(kr,kay2,3)
      kr=dsqrt(kay2)*dsin(beta)
      kth=-dsqrt(kay2)*dsin(pid2-beta)*dcos(az1)
      kph= dsqrt(kay2)*dsin(pid2-beta)*dsin(az1)
      IPHSPEED=1     ! 17 Mar 2006
      RETURN

C ----------------------------------------------------------------------!-------
C phase speed - 14 Mar 2006
      ENTRY PHSPEED(x)                                             ! 17 Mar 2006
      r=x
      call disper
! 16 Mar 2006      if (kz2ar.ge.0.d0) kzr=dsign(DSQRT(kz2ar),RCVRH)+kzbr     
c      if (kz2ar.lt.0.d0) kzr=0.0d0
c      print*,'phase kz2ar = ',kaz2r
      if (kz2ar.ge.0.d0)
     1 kz2r=2.d0*kzbr*dsign(DSQRT(kz2ar),RCVRH)+kz2ar+kzbr**2      ! 16 Mar 2006
      if (kz2ar.lt.0.d0) kz2r=kz2ar                                !  8 Apr 2011
      TEMPk2=kz2r+kth2+kph2                                        !  8 Apr 2011
      if (TEMPk2.gt.0.d0) PHSPEED=OW/DSQRT(TEMPk2)   ! 17 Mar 2006 !  8 Apr 2011
      if (TEMPk2.lt.0.d0) PHSPEED=-OW/DSQRT(-TEMPk2) ! 17 Mar 2006 !  8 Apr 2011
      if (TEMPk2.eq.0.d0) PHSPEED=0.d0                             !  8 Apr 2011
      RETURN

C ----------------------------------------------------------------------!-------
C Validity of ray theory                                           - 31 Mar 2006
      ENTRY IVRT(X,TLTP,TLTY,TLTL)
      CALL IDISPER
      WRITE(TLTP,'(2(1X,A8,F7.2,1X))') MODT
      WRITE(TLTL,'(10A4)') (TGP(I),I=1,10)
c      CALL ANNFIL('!P#hase !I#ntegral times 3root2 ('//            ! 6 Apr 2006
c     1'#RAD!)'         ,TLTY)
! 2 May 2006      CALL ANNFIL('!P#hase !I#ntegral (RAD) times 3root2!',TLTY)!10 Apr 2006 
      CALL ANNFIL('P#hase !I#ntegral (RAD) times 3root2!',TLTY)!10 Apr 2006 ! 2 May 2006
c      WRITE(TLTL,'(A40)') modrin

      r=earthr + dabs(rcvrh)
      th=pid2-llat
      ph=llon
      kr=dsin(beta)
      kth=-dsin(pid2-beta)*dcos(az1)
      kph= dsin(pid2-beta)*dsin(az1)
      call disper
      IVRT=0
      if(kay2.le.0.d0) return
c      call renorm(kr,kay2,3)
      kr=dsqrt(kay2)*dsin(beta)
      kth=-dsqrt(kay2)*dsin(pid2-beta)*dcos(az1)
      kph= dsqrt(kay2)*dsin(pid2-beta)*dsin(az1)
      kth2=kth*kth
      kph2=kph*kph
      IVRT=1
      xold = hb + earthr ! correction                              - 31 Mar 2006
c      print*,'subint = ',subint
      VRTT = subint                                                 ! 5 Apr 2006
      RETURN

C ----------------------------------------------------------------------!-------
C Validity of ray theory                                           - 31 Mar 2006
      ENTRY VRT(x)
      r=x
      call disper

      diff = x - xold
c      print*,'diff = ',diff

      if (kz2ar.ge.0.d0)
     1 VINT = (DSQRT(kz2ar) + kzbr)*3*DSQRT(2.d0)
      if (kz2ar.lt.0.d0)
     1 VINT = (-DSQRT(-kz2ar) + kzbr)*3*DSQRT(2.d0)

c      VINT = CABS(CSQRT(kz2a) + kzb)*3*DSQRT(2.d0)

      xold = x
      VRTT = VRTT + VINT*diff
c      print*,'VRTT = ',VRTT
      VRT = VRTT
      RETURN
C
C ----------------------------------------------------------------------!-------
      END FUNCTION PRFSUP                                          ! 27 Feb 2006
C ----------------------------------------------------------------------!-------

