C file S1.FOR
C ----------------------------------------------------------------------!-------
! 30 May 2013 - Corrected criterion for diagnostic printing
!  9 Mar 2012 - Replaced ' small rt' by ' X wave  '
!  9 Mar 2012 - Replaced ' usual rt' by ' O wave  '
!  9 Mar 2012 - Replaced ' root 3  ' by ' usual rt'
! 16 Sep 2011 - diagnostic printing only for GWWNL dispersion relation model
! 13 Apr 2011 - Modified putting a plus sign on raypaths when wanted
! 12 Apr 2011 - Added putting a plus sign on raypaths when wanted
! 14 Mar 2011 - Added diagnostic printing at transmitter
! 23 Feb 2011 - Changed criterion for diagnostic printing
! 21 Feb 2011 - Changed criterion for diagnostic printing
! 29 Nov 2010 - Altered the definition of diagnostic printing in sub. backup
! 22 Nov 2010 - Changed criterion for diagnostic printing
! 19 Nov 2010 - Changed criterion for diagnostic printing
! 19 Nov 2010 - Added errRell,errRelb,errRela to /rink/
!  5 Aug 2010 - Added kaysqi,kaysqbi,kaysqai to /rink/
! 29 Jul 2010 - Added /rink/
! 29 Jul 2010 - Added diagnostic printing
C 30 Apr 2009 - replaced REAL*8 with DOUBLE PRECISION
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 16 Mar 2006 - Added kzb to /RIN/
C 16 Mar 2006 - Changed kz2 to kz2a
C 14 Mar 2006 - Added kz2 to /RIN/
C  2 Mar 2006 - Moved CLEAR, ND2B, UCON to S5.for
C ----------------------------------------------------------------------!-------
Contains TRACE, PCROSS, RCROSS, HAMLTN, RKAM, RKAM1, BACKUP
C ----------------------------------------------------------------------!-------
C 17 Feb 2006 - Added testing to HAMLTN
C 16 Feb 2006 - Added integration of pressure to HAMLTN
C 2 Jun 2004 - fixed two bugs related to character string length
C ----------------------------------------------------------------------!-------
      SUBROUTINE TRACE
C ----------------------------------------------------------------------!-------
C     CALCULATES ONE RAYPATH FOR THE REQUESTED NUMBER OF HOPS
C     (CROSSINGS OR CLOSEST APPROACH TO RECEIVER SURFACE; A CLOSEST
C     APPROACH COUNTS AS TWO HOPS).
 
      implicit double precision (a-h,o-z)

 
      DOUBLE PRECISION NORYST
      PARAMETER (NORYST=0.0,MXPRGS=200)
C
C     COMMON DECK "SS" INSERTED HERE
      DOUBLE PRECISION MODSURF
      COMMON/SS/ MODSURF(4)
      COMMON/SS/U,PUR,PURR,PURTH,PURPH
      COMMON/SS/PUTH,PUPH,PUTHTH,PUPHPH,PUTHPH,USELECT,UTIME

C     COMMON DECK "RK" INSERTED HERE
C     DEFINE SIZE REQUIRED FOR RAY STATE SAVE ARRAY
      PARAMETER (LRKAMS=87+2*100,NXRKMS=12+LRKAMS,MXEQPT=21)
      PARAMETER (NRKSAV=NXRKMS+MXEQPT-1)
      COMMON /RK/ NEQS,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART

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 "TRAC" INSERTED HERE
      LOGICAL GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY
      COMMON /TRAC/ GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY,SMT,OSMT
      COMMON/TRAC/ROLD(20),DROLD(20),TOLD,ZDOT,D2Z,RAD,RAD1

C     COMMON DECK "RR" INSERTED HERE
      DOUBLE PRECISION MODREC
      COMMON/RR/ MODREC(4)
      COMMON/RR/F,PFR,PFRR,PFRTH,PFRPH
      COMMON/RR/PFTH,PFPH,PFTHTH,PFPHPH,PFTHPH,FSELECT,FTIME

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 "RKTIME" INSERTED HERE
      COMMON/CRKTIME/RKTIME

C     COMMON DECK "RINPLEX" INSERTED HERE
      DOUBLE PRECISION KAY2,KAY2I 
     & , kz2a,kz2ai,kzb,kzbi                                    ! 16 Mar 2006
      COMPLEX*16 PNP,POLAR,LPOLAR
      LOGICAL SPACE
      CHARACTER DISPM*6,modrin*64
      COMMON/RINPL/DISPM
      COMMON/RIN/MODRIN,RAYNAME(2,3),TYPE(3),SPACE
      COMMON/RIN/OMEGMIN,OMEGMAX,KAY2,KAY2I 
     & , kz2a,kz2ai,kzb,kzbi                ! 16 Mar 2006
      COMMON/RIN/PNP(10),POLAR,LPOLAR,SGN
      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
      dimension dd(0:4)
      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
     & , ks                                                       ! 20 Jul 2010

      COMMON R(20),T,STP,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))
     & ,(plusplot,W(78))                                           ! 12 Apr 2011
      LOGICAL PASTIME                                              ! 13 Apr 2011
C ----------------------------------------------------------------------!-------

C    THE FOLLOWING IS A LOCAL COMMON ONLY
C     COMMON DECK "TRLOCAL" INSERTED HERE
      COMMON/TRLOCAL/RSIGN,HOME,FDOT,GDOT,GOLD,GDOLD,UDOT,UOLD,UDOLD
      LOGICAL PCROSS,HOME,WASNT
      LOGICAL LAUNCH,APOGEE,PLTENB,EXTON
      EQUIVALENCE (R(2),TH),(R(3),PHI)
C
      CHARACTER*10 XCOND,DEFCND
      PARAMETER (DEFCND=' MAX HOPS')
C
      DOUBLE PRECISION PRGHST(3,MXPRGS)
      EQUIVALENCE   (APOG,PRIGEE,W(80))
      EQUIVALENCE (BOTABS,W(19))
C
c
      real *8 gtoo(12),ftoo(12),utoo(12)
      equivalence (gtoo(1),g),(ftoo(1),f),(utoo(1),u)
c
      LOGICAL PCNTRL                                               ! 29 Jul 2010
      PCNTRL=W(73).NE.0.d0                                         ! 29 Jul 2010
      PASTIME = PLUSPLOT.EQ.0.D0                                   ! 13 Apr 2011
      abc=setget(abc)
C     INCREMENT EVENT COUNT
c%      IRKTIME=IRKTIME+1
      rktime=rktime+1.d0
      PLTENB=APOG.EQ.0.d0
C     ENABLE EXTINCTION TEST IF LIMIT IS SET AND VARIABLE IS BEING INTE
      EXTON=(EXTINC.GT.0.d0) .AND. (W(58).NE.0.d0)
C     POINT INDEX TO ABSORTION VARIABLE
      NR=7
      IF(W(57).NE.0.d0) NR=8
      RMAX=HMAX+EARTHR
      RMIN=HMIN+EARTHR
      NHOP=HOP
      MAX=MAXSTP
      NSKIP=SKIP
C********* INITIALIZE PARAMETERS FOR INTEGRATION SUBROUTINE RKAM
      MODE=INTYP
      STEP=STEP1
      E1MAX=MAXERR
      E1MIN=MAXERR/ERATIO
      E2MAX=STPMAX
      E2MIN=STPMIN
      FACT=FACTR
      RSTART=1.d0
      CALL HAMLTN
      fdot=dot(ftoo)
c      FDOT=DOT(F)
C
C     CHECK FOR EQUALITY WITH RECEIVER HEIGHT WITHIN MACHINE PRECISION
      THERE=F.EQ.0.d0
C
      HOME=FDOT*F.GE.0.d0
C********* IHOP=0 TELLS PRINTR TO PRINT HEADING AND PUNCH A TRANSMITTER
C********* RAYSET AND NEWTRC=TRUE TELLS RAYPLT TO START A NEW RAY
      NEWTRC=.TRUE.
      IHOP=0
C********* RESET PERIGEE PLOT COUNTER
      NPRGS=0
C     USE CURRENT RELATIVE POSITION OF RAY TO PREDICT
C     SIGN OF NEXT INTERSECTION OF RAY WITH RECEIVER HEIGHT
C     (WILL SOON BE REVERSED).
c
      RSIGN=DSIGN(1.d0,F)
C     IF RAY IS LAUNCHED FROM RECEIVER HEIGHT, USE DIRECTION INSTEAD
      IF(THERE) RSIGN=DSIGN(1.d0,FDOT)
      if(there .and. fdot .eq. 0.d0) rsign=dsign(1.d0,d2z)
C  If the ray is also tangent to the receiver surface, use ray curvature
      gdot=dot(gtoo)
      udot=dot(utoo)
c      UDOT=DOT(U)
c      GDOT=DOT(G)
C     TEST CASE FOR RECEIVER ON THE GROUND
      IF(F.EQ.G) RSIGN=1.d0
C     IF RECEIVER IS ON TERRAIN THEN DIRECTION OF FIRST
C     RECEIVER CROSSING IS DOWNWARD
      CALL PRINTR('TXMTR    ',RAYSET)
      IF (PLTENB) CALL RAYPLT

      IF (W(73).eq.2.d0.or.W(73).eq.3.d0) 
     &  CALL PRINTR(' diagnose',NORYST)                            ! 14 Mar 2011
C ----------------------------------------------------------------------!-------

      NEWRAY=.TRUE.
      STH0=dSIN(TH)
      CTH0=dCOS(TH)
      PHI0=PHI
C
C     SET DEFAULT RAY EVENTS FOR 'PRINTR'
      XCOND=DEFCND
      XSET=NORYST
C
C********* LOOP ON NUMBER OF HOPS
10    IHOP=IHOP+1
C      REVERSE SIGN AT EACH CROSSING OF RECEIVER HEIGHT
      RSIGN=-RSIGN
      IF (IHOP.GT.NHOP) GO TO 100
      PENET=.FALSE.
C********* LOOP ON MAXIMUM NUMBER OF STEPS PER HOP
      DO 79 J=1,MAX
C     LAUNCH=TRUE ONLY WHEN TRANSMITTER ON GROUND WITH DOWNWARD
C     TRANSMISSION(POSSIBLE ONLY ON 1ST STEP).
      LAUNCH=G.EQ.0.d0 .AND. GDOT.LT.0.d0
C
C     SAVE CURRENT STATE VALUES FOR LATER COMPARISONS
C     RESTORE RAY STATE IF NEEDED
12    DO 13 L=1,6
      ROLD(L)=R(L)
13    DROLD(L)=DRDT(L)
      FDOLD=FDOT
      GDOLD=GDOT
      UDOLD=UDOT
      FOLD=F
      GOLD=G
      UOLD=U
      TOLD=T
C
C     PROCESS NEXT RAY POINT
      CALL RKAM
      fdot=dot(ftoo)
c      FDOT=DOT(F)
      WASNT=.NOT.HOME
      HOME=FDOT*F.GE.0.d0
C
C     CHOOSE STRATEGY A OR B
      IF(DROLD(1).GT.0.d0) GO TO 15
C
C     BEGIN STRATEGY A
C
C     LOOK FOR DOWNGOING CROSSING OF RECEIVER SURFACE
      IF(.NOT.LAUNCH.AND.F.LT.0.d0.AND.RSIGN.LT.0.d0) GO TO 50
C     CHECK FOR CASE OF CLOSEST APPROACH
      IF (HOME.AND.WASNT) GO TO 30
C     DETERMINE IF A GROUND CROSSING HAS OCCURED
      IF(PCROSS(gtoo,GDOT)) GO TO 20
C     LOOK FOR UPGOING CROSSING OF RECEIVER SURFACE
      IF(.NOT.LAUNCH.AND.F.GT.0.d0.AND.RSIGN.GT.0.d0) GO TO 50
C     CROSSED UPPER SURFACE?
      IF(PCROSS(utoo,UDOT)) GO TO 60
      GO TO 18
C
C     BEGIN STRATEGY B
C
C     LOOK FOR UPGOING CROSSING OF RECEIVER SURFACE
15    IF(.NOT.LAUNCH.AND.F.GT.0.d0.AND.RSIGN.GT.0.d0) GO TO 50
C     CHECK FOR CASE OF CLOSEST APPROACH
      IF (HOME.AND.WASNT) GO TO 30
C     CROSSED UPPER SURFACE?
      IF(PCROSS(utoo,UDOT)) GO TO 60
C     LOOK FOR DOWNGOING CROSSING OF RECEIVER SURFACE
      IF(.NOT.LAUNCH.AND.F.LT.0.d0.AND.RSIGN.LT.0.d0) GO TO 50
C     DETERMINE IF A GROUND CROSSING HAS OCCURED
      IF(PCROSS(gtoo,GDOT)) GO TO 20
C
C     RAY HAS NOT CROSSED ANY SURFACE OF DISCONTINUITY
C     OR RECEIVER SURFACE. PRINT OUT ANY SPECIAL CONDITIONS
C     OF THE RAY.
C
C     CHECK FOR PERIGEE CONDITION
18    IF (DROLD(1).GE.0.d0.OR.DRDT(1).LE.0.d0) GO TO 25
      CALL PRINTR(' PERIGEE ',NORYST)
      IF(PRIGEE.EQ.0.d0) GO TO 25
C
      IF(NPRGS.GE.MXPRGS) CALL STOPIT('PRG LIMIT')
      NPRGS=NPRGS+1
      CALL RMOVE(PRGHST(1,NPRGS),R,3)
C
25    APOGEE=DROLD(1).GT.0.d0.AND.DRDT(1).LT.0.d0
      PLTENB=APOG.EQ.0.d0.OR.APOGEE
C
      IF(APOGEE) CALL PRINTR(' APOGEE  ',NORYST)
      IF (DROLD(2)*DRDT(2).LT.0.d0) CALL PRINTR(' MAX LAT ',NORYST)
      IF (DROLD(3)*DRDT(3).LT.0.d0) CALL PRINTR(' MAX LONG',NORYST)
      DO 14 I=4,6
      IF (ROLD(I)*R(I).LT.0.d0) CALL PRINTR(' WAVE REV',NORYST)
   14 CONTINUE
      GO TO 75
C********* RAY WENT UNDERGROUND
C     USE 'FULL' REGULAR PARABOLIC FIT(ENTRY 'FIT')
20    IF(BOTABS.NE.0.d0) THEN
         CALL RCROSS(-1.d0,gtoo,GDOT,'XBOTM ABS',GROUND)
         XCOND=' '
      ELSE
         CALL RCROSS(-1.d0,gtoo,GDOT,'GBOTM REF',GROUND)
      ENDIF
C
      GO TO 75
C********* RAY WENT THROUGH UPPER SURFACE
C     USE 'FULL' REGULAR PARABOLIC FIT(ENTRY 'FIT')
60    CALL RCROSS(1.d0,utoo,UDOT,'ASURF REF',SURF)
      GO TO 75
C********* RAY MAY HAVE MADE A CLOSEST APPROACH
C     USE 'FULL' REGULAR PARABOLIC FIT
30    CALL FIT(ftoo,FOLD,FDOLD)
      IF(RAD.GE.0.d0) GO TO 501
C     ESTIMATE TIME OF CLOSEST APPROACH
      TP=T-ZDOT/D2Z
      CALL GRAZE(ftoo,RSIGN,TP)
      FDOT=ZDOT
      IF (THERE) GO TO 51
C
C     SET DRDT(1)=0 TO AVOID INCORRECT APOGEE OR PERIGEE PRINTOUT, A NEW
C     VALUE WILL BE COMPUTED BEFORE FURTHER ANALYSIS IS DONE
40    DRDT(1)=0.d0
      HPUNCH=R(1)-EARTHR
      CALL PRINTR('MMIN DIST',RAYSET)
      IF(PLTENB) CALL RAYPLT
      IF (IHOP.GE.NHOP) GO TO 100
      IHOP=IHOP+1
      RSIGN=-RSIGN
      CALL PRINTR ('MMIN DIST',RAYSET)
      GO TO 89
C********* RAY CROSSED RECEIVER HEIGHT
50    CALL FIT(ftoo,FOLD,FDOLD)
C     ESTIMATE GROUP TIME WHEN RAY CROSSES RECEIVER HEIGHT
C     (IN THE CORRECT DIRECTION).
501   TC=T-2.d0*F/(ZDOT+SIGN(RAD1,RSIGN))
      CALL BACKUP(ftoo,RSIGN,TC)
      FDOT=ZDOT
      IF(.NOT.THERE) GO TO 40
51    HPUNCH=R(1)-EARTHR
      CALL PRINTR('RRCVR    ',RAYSET)
      IF (PLTENB) CALL RAYPLT
      IF (GET(ftoo).NE.GET(gtoo)) GO TO 89
C
C     RECEIVER IS ON TERRAIN
      IF (IHOP.GE.NHOP) GO TO 100
      IHOP=IHOP+1
      RSIGN=-RSIGN
C********* GROUND REFLECT
      r777=REFLECT(gtoo)
      CALL HAMLTN
      FDOT=DOT(ftoo)
      RSTART=1.d0
      HPUNCH=R(1)-EARTHR
      CALL PRINTR('GBOTM REF',RAYSET)
      THERE=.TRUE.
      HPUNCH=R(1)-EARTHR
      CALL PRINTR ('RRCVR    ',RAYSET)
      GO TO 89
C*********
75    IF (PLTENB) CALL RAYPLT
C ----------------------------------------------------------------------!-------
      IF (PLTENB.and.T.GE.plusplot.AND..NOT.PASTIME) THEN 
         call rayplus                                              ! 12 Apr 2011
         PASTIME = .TRUE.                                          ! 13 Apr 2011
      ENDIF                                      
      PLTENB=APOG.EQ.0.d0
      IF(EXTON) THEN
         IF (R(NR).GE.EXTINC) XCOND='EEXTINC'
      ENDIF
      IF (R(1).GT.RMAX.AND.F.GT.0.d0.AND.FDOT.GT.0.d0) XCOND='UMAX HT'
      IF (R(1).LT.RMIN.AND.F.LT.0.d0.AND.FDOT.LT.0.d0) XCOND='DMIN HT'
      RANGE=EARTHR*dACOS(dCOS(TH)*CTH0+dSIN(TH)*STH0*dCOS(PHI-PHI0))
      IF(RGMAX.GT.0.d0.AND.RGMAX.LT.RANGE) XCOND='FMAX RANGE'
      IF(XCOND.NE.DEFCND) GO TO 90
C
      IF (MOD(J,NSKIP).EQ.0) CALL PRINTR('         ',NORYST)      !  2 Jun 2004

c      print *,'dispm =',dispm,'before first loop'
      IF (dispm.eq.'GWWNL ') THEN                                 ! 16 Sep 2011
c      print *,'dispm =',dispm,'in first loop'
!                                         corrected below         ! 30 May 2013
       IF ((W(73).eq.2.d0.or.W(73).eq.3.d0).and.ray.eq.0.d0.and.  ! 22 Nov 2010
     & (abs(errRell(n3)).gt.abs(errRell(nsmall)))   ! 19 Nov 2010 ! 21 Feb 2011
     &  .and.abs(errRell(nsmall)).NE.0.D0) then                   ! 23 Feb 2011
         CALL PRINTR(' X wave  ',NORYST)! 22 Nov 2010! 21 Feb 2011 ! 9 Mar 2012
         CALL PRINTR('         ',NORYST)                          ! 21 Feb 2011
       ENDIF                                                      ! 21 Feb 2011
      ENDIF                                                       ! 16 Sep 2011

c      print *,'dispm =',dispm,'before second loop'
      IF (dispm.eq.'GWWNL ') THEN                                 ! 16 Sep 2011
c      print *,'dispm =',dispm,'in second loop'
!                                         corrected below         ! 30 May 2013
       IF ((W(73).eq.2.d0.or.W(73).eq.3.d0).and.ray.eq.-1.d0.and. ! 22 Nov 2010
     & (abs(errRell(n3)).lt.abs(errRell(nsmall)))   ! 19 Nov 2010 ! 21 Feb 2011
     &  .and.abs(errRell(n3)).NE.0.D0) then                       ! 23 Feb 2011
         CALL PRINTR(' O wave  ',NORYST)!22 Nov 2010! 21 Feb 2011 !  9 Mar 2012
         CALL PRINTR('         ',NORYST)                          ! 21 Feb 2011
       ENDIF                                                      ! 21 Feb 2011
      ENDIF                                                       ! 16 Sep 2011
79    CONTINUE
C********* EXCEEDED MAXIMUM NUMBER OF STEPS
      HPUNCH=R(1)-EARTHR
      CALL PRINTR('SSTEP MAX',RAYSET)
      GO TO 100
C
C*********
89    HOME=.TRUE.
      GDOT=DOT(gtoo)
      GO TO 10
C********* RAY PENETRATED COMPUTATIONAL AREA BOUNDARY AND WAS HEADING OU
90    PENET=.TRUE.
      HPUNCH=R(1)-EARTHR
      XSET=RAYSET
C
C     NORMAL EXIT
100   IF(XCOND.NE.' ') CALL PRINTR(XCOND,XSET)
      IF(NPRGS.LE.0) RETURN
C
      CALL DASH
      NEWTRC=.TRUE.
      DO 110 I=1,NPRGS
      CALL RMOVE(R,PRGHST(1,I),3)
110   CALL RAYPLT
C
      CALL RESET('DASH')
C ----------------------------------------------------------------------!-------
      END SUBROUTINE TRACE
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      LOGICAL FUNCTION PCROSS(Z,ZDT)
C ----------------------------------------------------------------------!-------
C     RETURNS 'TRUE' IF PARABOLIC 'FIT2'(SEE 'FIT' ROUTINE)
C     INDICATES CROSSING OF SURFACE 'Z'.
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      DOUBLE PRECISION ZDT(3),z(12)
C     COMMON DECK "TRAC" INSERTED HERE
      LOGICAL GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY
      COMMON /TRAC/ GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY,SMT,OSMT
      COMMON/TRAC/ROLD(20),DROLD(20),TOLD,ZDOT,D2Z,RAD,RAD1
C
C     USE SPECIAL PARABOLIC FIT GIVING 'SMT' AND 'OSMT'
      ZDT(1)=DOT(Z)
      PCROSS=.TRUE.
      IF(Z(1).LT.0.d0) GO TO 20
      IF(ZDT(1).LE.0.d0.OR.ZDT(3).GE.0.d0) GO TO 15
      CALL FIT2(Z,ZDT(2),ZDT(3))
      IF(SMT.GT.Z(1).OR.OSMT.GT.ZDT(2)) GO TO 20
C
15    PCROSS=.FALSE.
20    RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION PCROSS
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE RCROSS(S,Z,ZDT,EVENT,QMODE)
C ----------------------------------------------------------------------!-------
C     FIND ESTIMATED CROSSING POINT OF SURFACE 'Z' THEN USE
C     ROUTINE 'BACKUP' TO GO THERE AND PERFORM A RAY REFLECTION
C     WITH ROUTINE 'REFLECT'.
C     EXTENDS RAY PROPAGATION TO INTERSECTION WITH REFLECTING SURFACE
C     Z AND CALLS 'REFLECT' TO OBTAIN VECTOR COMPONENTS.
C
 
      implicit double precision (a-h,o-z)
 
      CHARACTER EVENT*9 ! 2 Jun 2004
      DOUBLE PRECISION ZDT(3),z(12)
      LOGICAL QMODE
C     COMMON DECK "TRAC" INSERTED HERE
      LOGICAL GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY
      COMMON /TRAC/ GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY,SMT,OSMT
      COMMON/TRAC/ROLD(20),DROLD(20),TOLD,ZDOT,D2Z,RAD,RAD1
C     COMMON DECK "RR" INSERTED HERE
      DOUBLE PRECISION MODREC
      COMMON/RR/ MODREC(4)
      COMMON/RR/F,PFR,PFRR,PFRTH,PFRPH
      COMMON/RR/PFTH,PFPH,PFTHTH,PFPHPH,PFTHPH,FSELECT,FTIME
      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))
C     COMMON DECK "FLAG" INSERTED HERE
      LOGICAL NEWWR,NEWWP,NEWTRC,PENET
      COMMON /FLG/ NTYP,NEWWR,NEWWP,NEWTRC,PENET,LINES,IHOP,HPUNCH
      COMMON/FLGP/NSET
      COMMON R(20),T,STP,DRDT(20)
      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART
      LOGICAL HOME
C     COMMON DECK "TRLOCAL" INSERTED HERE
      COMMON/TRLOCAL/RSIGN,HOME,FDOT,GDOT,GOLD,GDOLD,UDOT,UOLD,UDOLD
C
      DOUBLE PRECISION ftoo(12)
      equivalence (ftoo(1),f)
      IF(QMODE) GO TO 60
      CALL FIT(Z,ZDT(2),ZDT(3))
      TC=T-2.d0*z(1)/(ZDOT-RAD1)
      CALL BACKUP(Z,-1.d0,TC)
C********* GROUND REFLECT
60    r777=REFLECT(Z)
      CALL HAMLTN
      ZDT(1)=DOT(Z)
      FDOT=DOT(ftoo)
      HOME=FDOT*F.GE.0.d0
      RSTART=1.d0
      HPUNCH=R(1)-EARTHR
      CALL PRINTR(EVENT,RAYSET)
      IF(F.NE.z(1)) RETURN
C     AVOID RECEIVER CROSSING AT THE TRANSMITTER IF RECEIVER ON TERRAIN
      RSIGN=S
      HOME=.TRUE.
C ----------------------------------------------------------------------!-------
      END SUBROUTINE RCROSS
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE HAMLTN
C ----------------------------------------------------------------------!-------
C     CALCULATES HAMILTONS EQUATIONS FOR RAY TRACING
C     AND OTHER QUANTITIES TO BE INTEGRATED
C     COMMON DECK "GG" INSERTED HERE
 
      implicit double precision (a-h,o-z)
 
      DOUBLE PRECISION MODG,modc
      common/rkams/xv(5),fv(4,20),yu(5,20),epm,alpha,mm
      common/cc/modc(4),cs,pcst,pcsr,pcsth,pcsph           !Rev 7/26/90
      COMMON/GG/MODG(4)
      COMMON/GG/G,PGR,PGRR,PGRTH,PGRPH
      COMMON/GG/PGTH,PGPH,PGTHTH,PGPHPH,PGTHPH,GSELECT,GTIME

C     COMMON DECK "PP" INSERTED HERE
      REAL(KIND=8) p,modp
      COMMON/PP/MODP(4),p,ppt,ppr,ppth,ppph

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 "RINREAL" INSERTED HERE
      LOGICAL SPACE
      DOUBLE PRECISION LPOLAR,LPOLRI,KPHK,KPHKI,KAY2,KAY2I 
     & , kz2a,kz2ai,kzb,kzbi           ! 16 Mar 2006
      CHARACTER DISPM*6,modrin*64
      COMMON/RINPL/DISPM
      COMMON /RIN/ MODRIN,RAYNAME(2,3),TYPE(3),SPACE
      COMMON/RIN/OMEGMIN,OMEGMAX,KAY2,KAY2I 
     & , kz2a,kz2ai,kzb,kzbi                ! 16 Mar 2006
     1,             H,HI,PHT,PHTI,PHR,PHRI,PHTH,PHTHI,PHPH,PHPHI
     2, PHOW,PHOWI,PHKR,PHKRI,PHKTH,PHKTI, PHKPH,PHKPI
     3 ,KPHK,KPHKI,POLAR,POLARI,LPOLAR,LPOLRI,SGN

      COMMON R(20),T,STP,DRDT(20)
      COMMON/WW/DUM(10),MAXW,W(1000)
      EQUIVALENCE (TH,R(2)),(PH,R(3)),(KR,R(4)),(KTH,R(5)),(KPH,R(6)),
     1 (DTHDT,DRDT(2)),(DPHDT,DRDT(3)),(DKRDT,DRDT(4)),(DKTHDT,DRDT(5)),
     2 (DKPHDT,DRDT(6)),(HMAX,W(26)),(OW,W(6))
     3 ,(W(554),H0) ! 16 Feb 2006
      DOUBLE PRECISION KR,KTH,KPH
      STH=dSIN(TH)
      CTH=dSIN(PID2-TH)
      RSTH=R(1)*STH
      RCTH=R(1)*CTH
      CALL DISPER
      if (w(999).eq.1.d0)write(*,'(4e18.11)')r(1)-w(1),th,t,stp
      if (w(999).eq.2.d0)write(*,'(5f10.3 )')(xv(i),i=1,5)
      if (w(999).eq.3.d0)write(*,'(2f10.4 )')r(1)-w(1),sqrt(cs) !Rev 
      DENPHC=1.d0/(PHOW*CREF)
      DRDT(1)=-PHKR*DENPHC
      DTHDT=-PHKTH*DENPHC/R(1)
      DPHDT=-PHKPH*DENPHC/RSTH
      DKRDT=PHR*DENPHC+KTH*DTHDT+KPH*STH*DPHDT
      DKTHDT=(PHTH*DENPHC-KTH*DRDT(1)+KPH*RCTH*DPHDT)/R(1)
      DKPHDT=(PHPH*DENPHC-KPH*STH*DRDT(1)-KPH*RCTH*DTHDT)/RSTH
      NR=6
C********* PHASE PATH
      IF (W(57).EQ.0.d0) GO TO 10
      NR=NR+1
      DRDT(NR)=-KPHK/PHOW/OW
C********* ABSORPTION
   10 IF (W(58).EQ.0.d0) GO TO 15
      NR=NR+1
      DRDT(NR)= 10.d0/ALN10*KPHK*KAY2I/(KR*KR+KTH*KTH+KPH*KPH)*DENPHC
C********* DOPPLER SHIFT
   15 IF (W(59).EQ.0.d0) GO TO 20
      NR=NR+1
      DRDT(NR)=-PHT*DENPHC/PIT2
C********* GEOMETRICAL PATH LENGTH
   20 IF (W(60).EQ.0.d0) GO TO 25
      NR=NR+1
      DRDT(NR)=dSQRT(PHKR*phkr+PHKTH*phkth+PHKPH*phkph)*dABS(DENPHC)
C********* TERRAIN FUNCTION AS COMPLETE INTEGRAL
25    IF(dABS(W(61))+dABS(W(62))+dABS(W(63))+dABS(W(64)).EQ.0.d0)
     1 GO TO 45
      CALL TOPOG
C
      IF(W(61).EQ.0.d0) GO TO 30
      NR=NR+1
      DRDT(NR)=PGR*DRDT(1)+PGTH*DTHDT+PGPH*DPHDT
C********* TERRAIN TIME DERIVATIVES AS COMPLETE INTEGRALS
   30 IF(W(62).EQ.0.d0) GO TO 35
      NR=NR+1
C     DERIVATIVE WITH RESPECT TO R
      DRDT(NR)=PGRR*DRDT(1)+PGRTH*DTHDT+PGRPH*DPHDT
   35 IF(W(63).EQ.0.d0) GO TO 40
      NR=NR+1
C     DERIVATIVE WITH RESPECT TO THETA
      DRDT(NR)=PGRTH*DRDT(1)+PGTHTH*DTHDT+PGTHPH*DPHDT
   40 IF(W(64).EQ.0.d0) GO TO 45
      NR=NR+1
C     DERIVATIVE WITH RESPECT TO PHI
c      print *, ' PGRPH = ' , PGRPH 
c      print *, ' PGTHPH = ' , PGTHPH 
c      print *, ' PGPHPH = ' , PGPHPH 
      DRDT(NR)=PGRPH*DRDT(1)+PGTHPH*DTHDT+PGPHPH*DPHDT
c      print *, 'drdt(nr) = ' , drdt(nr)
c      print *, ' nr = ' , nr
c      print *, 'r(nr) = ' , r(nr)
   45 IF(W(65).EQ.0.d0) GO TO 50 ! 16 Feb 2006
      NR=NR+1 ! 16 Feb 2006
C     DERIVATIVE of Pressure ! 16 Feb 2006
c      testing = ppr/p ! 17 Feb 2006
c      testing2 = testing + 1.d0/h0 ! 17 Feb 2006
c      print *,' ppr/p = ', testing,' ppr/p +1/H0 = ', testing2 ! 17 Feb 2006
c      print *,'ppr = ',ppr, 'p = ',p ! 17 Feb 2006
      DRDT(NR)=ppt/p+(1.d0/H0+ppr/p)*DRDT(1)+(ppth*DTHDT+ppph*DPHDT)/p !16Feb 06
c      print *, 'drdt(nr) = ' , drdt(nr)
c      print *, ' nr = ' , nr
c      print *, 'r(nr) = ' , r(nr)
C********* OTHER CALCULATIONS
50    CONTINUE
      RETURN
C ----------------------------------------------------------------------!-------
         END SUBROUTINE HAMLTN
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE RKAM
C ----------------------------------------------------------------------!-------
C     KEEPS TRACK OF INTERNAL INTEGRATION STEPS PERFORMED BY THE RKAM1
C     ROUTINE AND MAKES THEM AVAILABLE TO CALLING ROUTINES.
 
      implicit double precision (a-h,o-z)
 
      COMMON/RKAMS/XV(5),FV(4,20),YU(5,20),EPM,ALPHA,MM
C     COMMON DECK "RK" INSERTED HERE
C     DEFINE SIZE REQUIRED FOR RAY STATE SAVE ARRAY
      PARAMETER (LRKAMS=87+2*100,NXRKMS=12+LRKAMS,MXEQPT=21)
      PARAMETER (NRKSAV=NXRKMS+MXEQPT-1)
      COMMON /RK/ NEQS,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART
C     COMMON DECK "RKTIME" INSERTED HERE
      COMMON/CRKTIME/RKTIME
      COMMON Y(20),T,SPACE,DYDT(20)
      DOUBLE PRECISION RV(1)
      DOUBLE PRECISION YU
      EQUIVALENCE(RV,R)
C
      DOUBLE PRECISION SVBUF(NRKSAV)
      LOGICAL SAVED
C
      DATA SAVED/.FALSE./
C     PERFORM CLOCK ADVANCE(SEE 'GET' ROUTINE)
c%      IRKTIME=IRKTIME+1
      rktime=rktime+1.d0
C
      IF(RSTART.NE.0.d0 .OR. MODE.LE.2) GO TO 250
C
      NV=NV+1
      IF(NV.LE.4) GO TO 300
C
250   TOLD=T
C     PERFORM NUMERICAL INTEGRATION TO TIME 'T'
      CALL RKAM1
      IF(MODE.LE.2) GO TO 400
      IF(RSTART.NE.0.d0) GO TO 260
C     SEARCH FOR REQUESTED VALUE IN PIPELINE
C     TAKING SIGN OF 'STEP' INTO ACCOUNT.
      TUP=TOLD+.25d0*STEP
      DO 150 NV=1,3
      IF(SPACE.GT.0.d0.AND.TUP.LT.XV(NV)) GO TO 300
      IF(SPACE.LT.0.d0.AND.TUP.GT.XV(NV)) GO TO 300
150   CONTINUE
C     GOT IT ALREADY, RETURN
      GO TO 400
C
260   NV=1
C     RETRIEVE REQUESTED VALUES FROM PIPELINE
300   T=XV(NV)
      DO 350 I=1,NEQS
      Y(I)=YU(NV,I)
350   DYDT(I)=FV(NV,I)
C
C     STANDARD EXIT SEQUENCE
400   RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY RKSAVE(SVBUF)
C ----------------------------------------------------------------------!-------
C
      SVBUF(1)=NV
      SVBUF(2)=NEQS
      SVBUF(3)=MM
      SVBUF(4)=SPACE
      SVBUF(5)=MODE
      CALL RMOVE(SVBUF(6),E1MAX,6)
      CALL RMOVE(SVBUF(12),XV,LRKAMS)
      CALL RMOVE(SVBUF(NXRKMS),Y,MXEQPT)
C
      SAVED=.TRUE.
      RETURN
C ----------------------------------------------------------------------!-------
      ENTRY RKRSTR(SVBUF)
C ----------------------------------------------------------------------!-------
C
      IF(.NOT.SAVED) RETURN
C
      NV=SVBUF(1)
      NEQS=SVBUF(2)
      MM=SVBUF(3)
      SPACE=SVBUF(4)
      MODE=SVBUF(5)
      CALL RMOVE(E1MAX,SVBUF(6),6)
      CALL RMOVE(XV,SVBUF(12),LRKAMS)
      CALL RMOVE(Y,SVBUF(NXRKMS),MXEQPT)
C
      WRITE(3,*) 'RESTORING NV,NEQS,MM,ALPHA,T='
     1   ,NV,NEQS,MM,ALPHA,T
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE RKAM
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE RKAM1
C ----------------------------------------------------------------------!-------
C        NUMERICAL INTEGRATION OF DIFFERENTIAL EQUATIONS
C     THIS ROUTINE IS A MODIFICATION OF RKAMSUB, WHICH WAS WRITTEN
C     BY G.J. LASTMAN AND IS AVAILABLE THROUGH THE CDC CO-OP LIBRARY
C     AS 'D2 UTEX RKAMSUB'.
 
      implicit double precision (a-h,o-z)
 
      COMMON /RK/ NN,SPACE,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART
      COMMON Y(20),T,STEP,DYDT(20)
      COMMON/RKAMS/XV(5),FV(4,20),YU(5,20),EPM,ALPHA,MM
      DIMENSION DELY(4,20),BET(4)
      DOUBLE PRECISION YU
C
      IF (RSTART.EQ.0.d0) GO TO 1000
      LL=1
      MM=1
      IF (MODE.EQ.1) MM=4
      ALPHA=T
      EPM=0.d0
      BET(1)=0.5d0
      BET(2)=0.5d0
      BET(3)=1.0d0
      BET(4)=0.0d0
      STEP=SPACE
      R=19.d0/270.d0
      XV(MM)=T
      IF (E1MIN.LE.0.d0) E1MIN=E1MAX/50.d0
      IF (FACT.LE.0.d0) FACT=0.5d0
      CALL HAMLTN
      DO 320 I=1,NN
      FV(MM,I)=DYDT(I)
  320 YU(MM,I)=Y(I)
      RSTART=0.d0
      GO TO 1001
 1000 IF (MODE.NE.1) GO TO 2000
C
C          RUNGE-KUTTA
 1001 DO 1034 K=1,4
      DO 1350 I=1,NN
      DELY(K,I)=STEP*FV(MM,I)
      Z=YU(MM,I)
 1350 Y(I)=Z+BET(K)*DELY(K,I)
      T=BET(K)*STEP+XV(MM)
      CALL HAMLTN
      DO 1034 I=1,NN
 1034 FV(MM,I)=DYDT(I)
      DO 1039 I=1,NN
      DEL=(DELY(1,I)+2.d0*DELY(2,I)+2.d0*DELY(3,I)+DELY(4,I))/6.d0
 1039 YU(MM+1,I)=YU(MM,I)+DEL
      MM=MM+1
      XV(MM)=XV(MM-1)+STEP
      DO 1400 I=1,NN
 1400 Y(I)=YU(MM,I)
      T=XV(MM)
      CALL HAMLTN
      IF (MODE.EQ.1) GO TO 42
      DO 150 I=1,NN
  150 FV(MM,I)=DYDT(I)
      IF (MM.LE.3) GO TO 1001
C
C          ADAMS-MOULTON
 2000 DO 2048 I=1,NN
      DEL=STEP*(55.d0*FV(4,I)-59.d0*FV(3,I)+37.d0*FV(2,I)-9.d0*
     1  FV(1,I))/24.d0
      Y(I)=YU(4,I)+DEL
 2048 DELY(1,I)=Y(I)
      T=XV(4)+STEP
      CALL HAMLTN
      XV(5)=T
      DO 2051 I=1,NN
      DEL=STEP*(9.d0*DYDT(I)+19.d0*FV(4,I)-5.d0*FV(3,I)+FV(2,I))/24.d0
      YU(5,I)=YU(4,I)+DEL
 2051 Y(I)=YU(5,I)
      CALL HAMLTN
      IF (MODE.LE.2) GO TO 42
C
C          ERROR ANALYSIS
      SSE=0.d0
      DO 3033 I=1,NN
      EPSIL=R*dABS(Y(I)-DELY(1,I))
      IF (MODE.EQ.3.AND.Y(I).NE.0.d0) EPSIL=EPSIL/dABS(Y(I))
      IF (SSE.LT.EPSIL) SSE=EPSIL
 3033 CONTINUE
      IF (E1MAX.GT.SSE) GO TO 3035
      IF (dABS(STEP).LE.E2MIN) GO TO 42
      LL=1
      MM=1
      STEP=STEP*FACT
      GO TO 1001
 3035 IF (LL.LE.1.OR.SSE.GE.E1MIN.OR.E2MAX.LE.dABS(STEP)) GO TO 42
      LL=2
      MM=3
      XV(2)=XV(3)
      XV(3)=XV(5)
      DO 5363 I=1,NN
      FV(2,I)=FV(3,I)
      FV(3,I)=DYDT(I)
      YU(2,I)=YU(3,I)
 5363 YU(3,I)=YU(5,I)
      STEP=2.d0*STEP
      GO TO 1001
C
C          EXIT ROUTINE
   42 LL=2
      MM=4
      DO 12 K=1,3
      XV(K)=XV(K+1)
      DO 12 I=1,NN
      FV(K,I)=FV(K+1,I)
   12 YU(K,I)=YU(K+1,I)
      XV(4)=XV(5)
      DO 52 I=1,NN
      FV(4,I)=DYDT(I)
   52 YU(4,I)=YU(5,I)
      IF (MODE.LE.2) RETURN
      E=dABS(XV(4)-ALPHA)
      IF (E.LE.EPM+.25d0*STEP) GO TO 2000
      EPM=E
      RETURN
C ----------------------------------------------------------------------!-------
         END SUBROUTINE RKAM1
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE BACKUP(Z,RSIGN,TC)
C ----------------------------------------------------------------------!-------
C     MOVES THE RAY TO THE CLOSEST INTERSECTION WITH THE RECEIVER
C     OR TERRAIN SURFACE(VARIABLES 'FSELECT' OR 'GSELECT' IN LABELED
C     COMMONS /RR/ OR /GG/, RESPECTIVELY, TELL WHICH KIND OF SURFACE).
C
 
      implicit double precision (a-h,o-z)
 
      CHARACTER*9 NBAK,NGRAZ,NTRY
C
      PARAMETER (PRNZTL=0.5d-4,PRNDZTL=1.d-6)
C
      COMMON /RK/ N,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART
C     COMMON DECK "TRAC" INSERTED HERE
      LOGICAL GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY
      COMMON /TRAC/ GROUND,SURF,PERIGE,THERE,MINDIS,NEWRAY,SMT,OSMT
      COMMON/TRAC/ROLD(20),DROLD(20),TOLD,ZDOT,D2Z,RAD,RAD1
      COMMON R(20),T,STP,DRDT(20)
      COMMON/WW/DUM(10),MAXW,W(1000)
      EQUIVALENCE (EARTHR,W(1)),(INTYP,W(41)),(STEP1,W(44))
      DOUBLE PRECISION KR,kparlel
      EQUIVALENCE (KR,R(4)),(DKRDT,DRDT(4))
      DOUBLE PRECISION INTYP
      LOGICAL PCNTRL
      DATA NBAK,NGRAZ/' BACK UP0', ' GRAZE 1'/
C
      NTRY=NBAK
      GO TO 100
C
C ----------------------------------------------------------------------!-------
      ENTRY GRAZE(Z,RSIGN,TC)
C ----------------------------------------------------------------------!-------
      NTRY=NGRAZ
C
C     DEFINE BASE STEP SIZE
100   STPB=STP
C
C     DEFINE BACKUP LOCATION TOLERANCES BASED ON INTEGRATION MODE
C
      IF(MODE.LT.3) THEN
         TOL=STEP1
      ELSEIF(MODE.EQ.3) THEN
         TOL=dABS(E1MAX*STPB)
      ELSE
         TOL=E1MAX
      ENDIF
C     USE THIS OR MINIMUM STEP SIZE WHICHEVER IS LESS
      TOL=DMIN1(TOL,FACT*E2MIN)
C
      PCNTRL=W(73).EQ.1.d0                                         ! 29 Nov 2010
      THERE=.TRUE.
C      STEP TO ESTIMATED CROSSING AT TC.
C********* DIAGNOSTIC PRINTOUT
      IF(PCNTRL) CALL PRINTR (NTRY,0.d0)
      STEP=TC-T
      STEP=DSIGN(DMIN1(DABS(STPB),DABS(STEP)),STEP)
      MODE=1
      RSTART=1.d0
      TOLD=T
      CALL RMOVE(DROLD,DRDT,3)
      ZDOLD=DOT(Z)
c diagnostic      print *,'after dot'
      CALL RKAM
c diagnostic      print *,'after rkam'
      RSTART=1.d0
      ZDOT=DOT(Z)
c diagnostic      print *,'after dot2'
      IF(NTRY.EQ.NGRAZ .OR. RSIGN*ZDOT.LE.0.d0) GO TO 12
C
C********* FIND NEAREST INTERSECTION OF RAY WITH THE HEIGHT HS
      DO 10 I=1,10
      STEP=-Z/ZDOT
      STEP=dSIGN(DMIN1(dABS(STPB),DABS(STEP)),STEP)
      IF (DABS(Z).LT.PRNZTL .AND. DABS(STEP).LT.TOL) GO TO 60
C********* DIAGNOSTIC PRINTOUT
      IF(PCNTRL) CALL PRINTR(' BACK UP1',0.d0)
      MODE=1
      RSTART=1.d0
      TOLD=T
      CALL RMOVE(DROLD,DRDT,3)
      ZDOLD=ZDOT
      CALL RKAM
      ZDOT=DOT(Z)
   10 RSTART=1.d0
C
C********* FIND NEAREST CLOSEST APPROACH OF RAY TO THE HEIGHT HS
12    THERE=.FALSE.
      DO 20 I=1,10
C     DO 'LOCAL' PARABOLIC FIT
      CALL FIT3(Z,ZOLD,ZDOLD)
      STEP=-ZDOT/D2Z
      STEP=DSIGN(DMIN1(DABS(STPB),DABS(STEP)),STEP)
      IF (DABS(ZDOT).LE.PRNDZTL .AND. DABS(STEP).LT.TOL) GO TO 60
C********* DIAGNOSTIC PRINTOUT
      IF(PCNTRL) CALL PRINTR (' GRAZE 1 ',0.d0)
      MODE=1
      RSTART=1.d0
      TOLD=T
      CALL RMOVE(DROLD,DRDT,3)
      ZOLD=Z
      ZDOLD=ZDOT
      CALL RKAM
      RSTART=1.d0
      IF (D2Z*Z.LT.0.d0) GO TO 30
      IF(KPARLEL(Z).EQ.0.d0) GO TO 60
   20 CONTINUE
      WRITE(3,350)
350   FORMAT(' ****** COULDN''T FIND CLOSEST APPROACH IN 10 STEPS')
      GO TO 60
C
   30 CONTINUE
C********* DIAGNOSTIC PRINTOUT
      IF(PCNTRL) CALL PRINTR (' BACK UP2',0.d0)
      MODE=1
C********* ESTIMATE DISTANCE TO NEAREST INTERSECTION GOING THE RIGHT
C********* DIRECTION OF RAY WITH HEIGHT HS
      CALL FIT3(Z,ZOLD,ZDOLD)
      STEP=(-ZDOT+RSIGN*RAD1)/D2Z
      RSTART=1.d0
      CALL RKAM
      RSTART=1.d0
C********* FIND NEAREST INTERSECTION OF RAY WITH HEIGHT HS
      DO 40 I=1,10
      ZDOT=DOT(Z)
      STEP=-Z/ZDOT
      STEP=DSIGN(DMIN1(DABS(STPB),DABS(STEP)),STEP)
      IF (DABS(Z).LT.PRNZTL .AND. DABS(STEP).LT.TOL) GO TO 50
C********* DIAGNOSTIC PRINTOUT
      IF(PCNTRL) CALL PRINTR (' BACK UP3',0.d0)
      MODE=1
      RSTART=1.d0
      CALL RKAM
   40 RSTART=1.d0
   50 THERE=.TRUE.
C********* RESET STANDARD MODE AND INTEGRATION TYPE
   60 MODE=INTYP
      STEP=STPB
      RETURN
C ----------------------------------------------------------------------!-------
         END SUBROUTINE BACKUP
C ----------------------------------------------------------------------!-------
 
