C file S2.FOR
C ----------------------------------------------------------------------!-------
!  9 Aug 2016 - The formulas used to change the normal component of the 
!               wave-normal vector at reflection that include Coriolis terms 
!               work only when the reflecting surface is horizontal.
!  9 Aug 2016 - Altered slightly the reflection formulas.
!  5 Aug 2016 - Added common block /rink3/ to help with reflection of wave
!  9 Mar 2012 - Replaced ' small rt' by ' X wave  '
!  9 Mar 2012 - Replaced ' usual rt' by ' O wave  '
!  2 Mar 2012 - Replaced ' root 3  ' by ' usual rt'
! 25 Aug 2011 - Changes recommended by Larry Baker
! 23 Aug 2011 - Added 1 column to on-screen printing of elev. angle 
! 23 Aug 2011 - Corrected printing in subroutine PRINTR
! 15 Mar 2011 - Added unconditional diagnostic printing for w73=4
! 14 Mar 2011 - Added diagnostic printing at transmitter
! 22 Feb 2011 - Modified diagnostic printing
! 21 Feb 2011 - Modified diagnostic printing
! 29 Nov 2010 - Modified diagnostic printing
! 22 Nov 2010 - Modified diagnostic printing
! 19 Nov 2010 - Added printing of errRell,errRelb,errRela
! 19 Nov 2010 - Added errRell,errRelb,errRela to /rink/
!  5 Aug 2010 - Added kaysqi,kaysqbi,kaysqai to /rink/
!  4 Aug 2010 - Modified diagnostic printing
! 29 Jul 2010 - Modified diagnostic printing
! 28 Jul 2010 - Added alpha,dd to /rink/
! 28 Jul 2010 - Modified diagnostic printing
! 27 Jul 2010 - Added printing of Hrel,errRel 
! 27 Jul 2010 - Added Hrel,errRel to /rink/
! 20 Jul 2010 - Added printing of ks and kay2
! 15 Jul 2010 - Added printing of n3, n4, nbig, nsmall
! 15 Jul 2010 - Added n3,n4,nbig,nsmall to common block /rink/
! 28 Jun 2010 - Added printing of relerr,relerra,relerrb 
! 28 Jun 2010 - Added relerr,relerra,relerrb to common block /rink/
! 15 Jun 2010 - Print more quartic root of k
! 14 Jun 2010 - Print more quartic root of k squared
! 11 Jun 2010 - Print quartic root of k squared
! 28 May 2010 - Put error into drdt(19)
C 30 Apr 2009 - replaced REAL*8 with DOUBLE PRECISION
C 11 Mar 2009 - Moved CONBLK to separate file CONBLK2.FOR
C 27 Aug 2008 - Commented out "call mach"
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 15 May 2007 - Added conditional printing of sound speed to PRINTR
C 20 Apr 2007 - Added log density to transmitter and receiver raysets
C 27 Mar 2007 - Added conditional printing of density to PRINTR
C  5 Mar 2007 - Added molecular fractions to common block /MM/
C 30 Apr 2006 - Corrected printing error in SUBROUTINE PRINTR
C 24 Mar 2006 - Moved GET, WCHANGE, and RENORM to S5
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 GET1, ITEST, ITOC, SETTRC, SET2, ISET2 to S5.for
C ----------------------------------------------------------------------!-------
Contains REFLECT, FIT, PRINTR
C ----------------------------------------------------------------------!-------
C 21 Feb 2006 - Correction
C 17 Feb 2006 - Added conditional printing of pressure to PRINTR
C 16 Feb 2006 - Added column heading for pressure
C 21 December 2005 - Modified printing of variables to coincide with EIGEN
C  2 Jun 2004 - SUBROUTINE CONBLK - Hollerith constants to character strings
C 10 Jun 2003 - line 341 (FUNCTION GET) added ,PK(10)
C 10 Jun 2003 - added line 634 (SUBROUTINE CONBLK) fixes data statement for CUEF
C 10 Jun 2003 - lines 668-674 (SUBROUTINE CONBLK) fixed some Hollerith constants
C 10 Jun 2003 - lines 678-693 (SUBROUTINE CONBLK) Holl. const. to char strings
C 10 Jun 2003 - line 837 (SUBROUTINE SET2) ENTRY ISET2 to SUB ISET2 & rewrote
C ----------------------------------------------------------------------!-------
      FUNCTION REFLECT(Z)
C ----------------------------------------------------------------------!-------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
C     COMPUTES NORMAL AND PARALLEL COMPONENTS OF THE K-VECTOR AT
C     REFLECTION POINTS TO A SURFACE. WIND EFFECTS ARE INCLUDED.
C
      DOUBLE PRECISION Z(12),KPARLEL,KNORM

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

      integer*2 nrow,ncol

      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 "RKAM" INSERTED HERE
      DOUBLE PRECISION KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)

      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 "FNDER" INSERTED HERE
      COMMON/FNDER/GY(20),NZ,NPZR,NPZRR,NPZRTH,NPZRPH,NPZTH
      COMMON/FNDER/NPZPH,NPZTHTH,NPZPHPH,NPZTHPH,NSELECT,NTIME
C
      DOUBLE PRECISION NR,NTH,NPH,KRNR,KPARR,KPARTH,KPARPH
      real *8 anr(3)
      EQUIVALENCE (anr(1),nr),(anr(2),nth),(anr(3),nph)
C
      IENTRY=1
      REFLECT=0.d0
      GO TO 5
C
C ----------------------------------------------------------------------!-------
      ENTRY KPARLEL(Z)
C ----------------------------------------------------------------------!-------
      KPARLEL=0.d0
      IENTRY=2
      GO TO 5
C
C ----------------------------------------------------------------------!-------
      ENTRY KNORM(Z)
C ----------------------------------------------------------------------!-------
      KNORM=0.d0
      IENTRY=3
C
5     zholdr=get(z)
      Z(1)=zholdr
      NR=Z(NPZR)
      NTH=Z(NPZTH)/R
      NPH=Z(NPZPH)/(R*dSIN(TH))
C
      r999=renorm(anr,1.d0,3)
C      CALL RENORM(NR,1.d0,3)
C     COMPUTE THE NORMAL COMPONENT OF K-VECTOR TO SURFACE
      KRNR=KR*NR+KTH*NTH+KPH*NPH
      IF(IENTRY.NE.3) GO TO 8
C     IF ENTRY 3 THEN WE ARE DONE
      REFLECT=KRNR
      RETURN
C
C     COMPUTE THE PARALLEL VECTOR COMPONENT
8     KPARR=KR-KRNR*NR
      KPARTH=KTH-KRNR*NTH
      KPARPH=KPH-KRNR*NPH
C
      IF(IENTRY.NE.2) GO TO 10
      REFLECT=dABS(KPARR)+dABS(KPARTH)+dABS(KPARPH)
      RETURN
C
10    OWIPAR=OW-KPARR*VR-KPARTH*VTH-KPARPH*VPH
      VNR=VR*NR+VTH*NTH+VPH*NPH
C
      print *, ' incident kr  = ', kr
      print *, ' incident kth = ', kth
      print *, ' incident kph = ', kph

      sq=kr**2+kth**2+kph**2                       ! Corrected a typo 5 Aug 2016
!      print *, ' sq  = ', sq

!      print *, ' nr  = ', nr
!      print *, ' nth = ', nth
!      print *, ' nph = ', nph
!      print *, ' vnr = ', vnr

!      FCTR=2.d0*(KRNR+OWIPAR*VNR/(CS-VNR*VNR))
! There should never be a component of the wind normal to the surface 
! at the surface.
      fctr=2.d0*krnr
!      print *, ' krnr = ', krnr

! This formula for fctr is correct only for a horizontal surface   !  5 Aug 2016
!      print *, ' ows = ', ows
!      print *, ' azz = ', azz
!      print *, ' bzz = ', bzz
!      print *, ' czz = ', czz
!      print *, ' czz0a = ', czz0a
!      print *, ' czz0b = ', czz0b
!      print *, ' czzf0 = ', czzf0
!      print *, ' czzfth = ', czzfth
!      print *, ' czzfph = ', czzfph
!      print *, ' czzfthth = ', czzfthth
!      print *, ' czzfphph = ', czzfphph
!      print *, ' czzfthph = ', czzfthph                            !  5 Aug 2016

      rad=azz*czz+bzz*bzz
!      print *, ' rad = ', rad
      root=0.d0
      if(rad.gt.0.d0) root=sqrt(rad)                               !  9 Aug 2016
      compareplus=(bzz+root)/azz
      compareminus=(bzz-root)/azz
      if(rad.gt.0.d0) fctr=2.d0*root/azz*sign(1.d0,kr-bzz/azz)     !  9 Aug 2016
      KR=KR-FCTR*NR
      KTH=KTH-FCTR*NTH
      KPH=KPH-FCTR*NPH

      sq=kr**2+kth**2+kph**2                       ! Corrected a typo 5 Aug 2016
!      print *, ' sq  = ', sq

      print *, ' compareplus  = ', compareplus
      print *, ' compareminus  = ', compareminus
      print *, ' reflected kr  = ', kr
      print *, ' reflected kth = ', kth
      print *, ' reflected kph = ', kph

!      drdt(20)=23.d0                  ! for test printing in gwwnl 27 May 2010

C ----------------------------------------------------------------------!-------
      END FUNCTION REFLECT
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE FIT(Z,ZOLD,ZDOLD)
C ----------------------------------------------------------------------!-------
C     COMPUTES THREE TYPES OF PARABOLIC FITS TO RAY PATH RELATIVE
C     TO TERRAIN.
C
 
      implicit double precision (a-h,o-z)
 
      DOUBLE PRECISION Z(12)
      DOUBLE PRECISION KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)
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 "FNDER" INSERTED HERE
      COMMON/FNDER/GY(20),NZ,NPZR,NPZRR,NPZRTH,NPZRPH,NPZTH
      COMMON/FNDER/NPZPH,NPZTHTH,NPZPHPH,NPZTHPH,NSELECT,NTIME
C
      DOUBLE PRECISION D2(3)
C
C     USE FIT OF APPENDIX 'J' OF REPORT 'WPL-103' WHICH USES
C     WEIGHTED ESTIMATE OF 1ST DERIVATIVE.
C
      IENTRY=1
      GO TO 5
C
C     USE MODIFIED FIT REQUIRING HEIGHTS OF PARABOLA VERTICES FROM
C     CURRENT AND PREVIOUS RAY POINTS.
C ----------------------------------------------------------------------!-------
      ENTRY FIT2(Z,ZOLD,ZDOLD)
C ----------------------------------------------------------------------!-------
      IENTRY=2
      GO TO 5
C
C     USE FIT OF APPENDIX U(LOCAL VALUE OF 1ST DERIVATIVE)
C ----------------------------------------------------------------------!-------
      ENTRY FIT3(Z,ZOLD,ZDOLD)
C ----------------------------------------------------------------------!-------
      IENTRY=3
C
5     ZDOT=DOT(Z)
C
      DTI=1.d0/(TPULSE-TOLD)
      DO 10 I=1,3
10    D2(I)=(DRDT(I)-DROLD(I))*DTI
C
      D2Z=Z(NPZR)*D2(1)+Z(NPZTH)*D2(2)+Z(NPZPH)*D2(3)
     1        +Z(NPZRR)*DRDT(1)*DRDT(1)
     1        +Z(NPZTHTH)*DRDT(2)*DRDT(2)
     1        +Z(NPZPHPH)*DRDT(3)*DRDT(3)
     1   +2.d0*(Z(NPZRTH)*DRDT(1)*DRDT(2)
     1        +Z(NPZRPH)*DRDT(1)*DRDT(3)
     1        +Z(NPZTHPH)*DRDT(2)*DRDT(3) )
C
C     THE STATEMENTS FROM HERE TO 'END FIT' IMPLEMENT THE
C     PARABOLIC FITS IN EQUATIONS J.1 AND U.3 OF THE TEXT.
C
      IF(IENTRY.NE.2) GO TO 30
      SMT=0.d0
      IF (D2Z.NE.0.d0) SMT=0.5d0*ZDOT*ZDOT/D2Z
C     USE FIT U.3 AT THE PREVIOUS POINT OF RAY PATH
      OSMT=0.d0
      IF (D2Z.NE.0.d0) OSMT=0.5d0*ZDOLD*ZDOLD/D2Z
      GO TO 2000
C
C     IMPLEMENTATION OF FIT FOR EQUATION J.1
30    IF(IENTRY.EQ.3) GO TO 1000
      ZDOTM=.5d0*(ZDOT+ZDOLD)
C     IMPLEMENT TESTS OF EQUATIONS J.2 AND J.3
      IF(dABS(ZDOTM).LE.5.d-2*dABS(ZDOT)) GO TO 1000
      FCTR=(Z(NZ)-ZOLD)*DTI/ZDOTM
      D2Z=FCTR*D2Z
      ZDOT=FCTR*ZDOT
C     END FIT
C
C     COMMON CODE FOR FIT AND FIT3
1000  RAD=ZDOT*ZDOT-2.d0*Z(NZ)*D2Z
      RAD1=DSQRT(DMAX1(RAD,0.d0))
C
2000  CONTINUE
C
C ----------------------------------------------------------------------!-------
      END SUBROUTINE FIT
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE CONBLK_test                                       ! 11 Mar 2009
C ----------------------------------------------------------------------!-------
C     DATA INITIALIZATION AND FILE OPENING SERVICE ROUTINE
C
C     COMMON DECK "RMACH" INSERTED HERE
 
      implicit double precision (a-h,o-z)
 
      COMMON/CRMACH/rmach(9)
C     COMMON DECK "HDR" INSERTED HERE
      CHARACTER*10 INITID*80,DAT,TOD,dinpname*12   !Rev 7/13/90
      COMMON/HDR/SEC,stsec                         !Rev 7/13/90
      COMMON/HDRC/INITID,DAT,TOD
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     COMMON DECK "CC" INSERTED HERE
      DOUBLE PRECISION MODC
      COMMON/CC/MODC(4),CS,PCST,PCSR,PCSTH,PCSPH
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 "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 "TT" INSERTED HERE
      DOUBLE PRECISION MODT
      COMMON/TT/MODT(4), T,PTT,PTR,PTTH,PTPH
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
C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN
C     COMMON DECK "CUCON" INSERTED HERE
      COMMON/UCONV/CNVV(5,4)
      CHARACTER PCV*3,CNVC*2
      COMMON/UCONC/PCV(5),CNVC(5,4)
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 "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
      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 "B9" INSERTED HERE
      INTEGER GMX,GNTBL,GITBL,GFRMTBL,IDSG(20)
      COMMON/B9/GMX,GNTBL(10),GITBL(10),GFRMTBL(10),GGP(113)
      EQUIVALENCE (GGP,IDSG),(ANG,GGP(11))
C     COMMON DECK "B2" INSERTED HERE
      INTEGER DUMX,DUNTBL,DUITBL,DUFRMTB,IDSDU(20)
      COMMON/B2/DUMX,DUNTBL(10),DUITBL(10),DUFRMTB(10),DUGP(10)
      EQUIVALENCE (DUGP,IDSDU)
C     COMMON DECK "B4" INSERTED HERE
      INTEGER DCMX,DCNTBL,DCITBL,DCFRMTB,IDSDC(20)
      COMMON/B4/DCMX,DCNTBL(10),DCITBL(10),DCFRMTB(10),DCGP(10)
      EQUIVALENCE (DCGP,IDSDC)
C     COMMON DECK "B6" INSERTED HERE
      INTEGER DTMX,DTNTBL,DTITBL,DTFRMTB,IDSDT(20)
      COMMON/B6/DTMX,DTNTBL(10),DTITBL(10),DTFRMTB(10),DTGP(10)
      EQUIVALENCE (DTGP,IDSDT)
C     COMMON DECK "B1" INSERTED HERE
      INTEGER UMX,UNTBL,UITBL,UFRMTBL,IDSU(20)
      COMMON/B1/UMX,UNTBL(10),UITBL(10),UFRMTBL(10),UGP(10)
      EQUIVALENCE (UGP,IDSU)
C     COMMON DECK "B3" INSERTED HERE
      INTEGER CMX,CNTBL,CITBL,CFRMTBL,IDSC(20)
      COMMON/B3/CMX,CNTBL(10),CITBL(10),CFRMTBL(10),CGP(512)
      EQUIVALENCE (CGP,IDSC),(ANC,CGP(11))
C     COMMON DECK "B5" INSERTED HERE
      INTEGER TMX,TNTBL,TITBL,TFRMTBL,IDST(20)
      COMMON/B5/TMX,TNTBL(10),TITBL(10),TFRMTBL(10),TGP(262)
      EQUIVALENCE (TGP,IDST),(ANT,TGP(11))
C     COMMON DECK "B7" INSERTED HERE
      INTEGER MMX,MNTBL,MITBL,MFRMTBL,IDSM(20)
      DOUBLE PRECISION MGP
      COMMON/B7/MMX,MNTBL(10),MITBL(10),MFRMTBL(10),MGP(10)
      EQUIVALENCE (MGP,IDSM)

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

C     COMMON DECK "MM" INSERTED HERE
      DOUBLE PRECISION M,MODM
      COMMON/MM/MODM(4),M,PMT,PMR,PMTH,PMPH
     1 ,XO2,XN2,XO,XN,XH2O,XO3 ! 5 Mar 2007

C     COMMON DECK "PP" INSERTED HERE
      DOUBLE PRECISION MODP
      COMMON/PP/MODP(4),P,PPT,PPR,PPTH,PPPH

C     COMMON DECK "AA" INSERTED HERE
      DOUBLE PRECISION MODA
      DOUBLE PRECISION MU,MUPT,MUPR,MUPTH,MUPPH
      DOUBLE PRECISION KAP,KAPPT,KAPPR,KAPPTH,KAPPPH
      COMMON/AA/MODA(4),MU,MUPT,MUPR,MUPTH,MUPPH
      COMMON/AA/KAP,KAPPT,KAPPR,KAPPTH,KAPPPH

C     COMMON DECK "CB11" INSERTED HERE
      INTEGER SMX,SNTBL,SITBL,SFRMTBL,IDSS(10)
      COMMON/B11/SMX,SNTBL(10),SITBL(10),SFRMTBL(10),SGP(11)
      EQUIVALENCE (SGP,IDSS),(ANS,SGP(11))

C     COMMON DECK "CB12" INSERTED HERE
      INTEGER DSMX,DSNTBL,DSITBL,DSFRMTB,IDSDS(10)
      COMMON/B12/DSMX,DSNTBL(10),DSITBL(10),DSFRMTB(10),DSGP(11)
      EQUIVALENCE (DSGP,IDSDS),(ANDS,DSGP(11))

C     COMMON DECK "LL" INSERTED HERE
      DOUBLE PRECISION MODL
      COMMON/LL/MODL(4),APH,APHPT,APHPR,APHPTH,APHPPH

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

C     COMMON DECK "B10" INSERTED HERE
      INTEGER DGMX,DGNTBL,DGITBL,DGFRMTB,IDSDG(20)
      COMMON/B10/DGMX,DGNTBL(10),DGITBL(10),DGFRMTB(10),DGGP(10)
      EQUIVALENCE (DGGP,IDSDG)

C     COMMON DECK "B8" INSERTED HERE
      INTEGER RMX,RNTBL,RITBL,RFRMTBL,IDSR(20)
      COMMON/B8/RMX,RNTBL(10),RITBL(10),RFRMTBL(10),RGP(10)
      EQUIVALENCE (RGP,IDSR)

C     COMMON DECK "CB17" INSERTED HERE
      INTEGER VMX,VNTBL,VITBL,VFRMTBL,IDSV(20)
      COMMON/B17/VMX,VNTBL(10),VITBL(10),VFRMTBL(10),VGP(53)
      EQUIVALENCE (VGP,IDSV),(ANV,VGP(11))

C     COMMON DECK "CB18" INSERTED HERE
      INTEGER DVMX,DVNTBL,DVITBL,DVFRMTB,IDSDV(20)
      COMMON/B18/DVMX,DVNTBL(10),DVITBL(10),DVFRMTB(10),DVGP(11)
      EQUIVALENCE (DVGP,IDSDV),(ANDV,DVGP(11))

C     COMMON DECK "CB19" INSERTED HERE
      INTEGER PRMX,PRNTBL,PRITBL,PRFRMTB,IDSPR(20)
      COMMON/B19/PRMX,PRNTBL(10),PRITBL(10),PRFRMTB(10),PRGP(11)
      EQUIVALENCE (PRGP,IDSPR),(ANP,PRGP(11))

C     COMMON DECK "CB20" INSERTED HERE
      INTEGER DPMX,DPNTBL,DPITBL,DPFRMTB,IDSDP(20)
      COMMON/B20/DPMX,DPNTBL(10),DPITBL(10),DPFRMTB(10),DPGP(11)
      EQUIVALENCE (DPGP,IDSDP),(ANDP,DPGP(11))
      COMMON/KNKN/KNBP,KNVC,KNDT
      COMMON/DDLIM/MXIX,MXIY,MNIX,MNIY
C
      DOUBLE PRECISION KVECT(22)
      DOUBLE PRECISION VSET(20)
      EQUIVALENCE(KVECT,KAY2),(VSET,V)
C
      DOUBLE PRECISION CUEF(3)
      INTEGER NUYIND(5)
C
      CHARACTER PFV(5)*3,CQVC(5,4)*2
C this sets the pconst common cref, rgas and gamma
C      DATA CUEF/1.d0,1.4d0,8.31436d-3/
	DATA CUEF/1.d0,8.31436d-3,1.4d0/ ! 10 Jun 2003
      DATA NUYIND/8,4,0,10,11/
      DATA PFV/'AN','LN','FQ','AM',' '/
      DATA CQVC/'RD' ,'KM' ,'RD' ,'DB',' '
     1        ,'DG' ,' M' ,'HZ' ,'NP',' '
     2        ,'KM' ,'FT' ,'KH' ,2*' '
     3        ,'  ' ,'NM' ,'DG' ,2*' '/
C
      CALL SET2(KVECT,0.d0,22)
      dummy=' ' ! 2Jun 2004
      CALL SET2(RAYNAME,dummy, 6)
      CALL SET2(VSET,0.d0,20)
      CALL SET2(OMEGMIN,0.d0,2)
      MXIX=-1000
      MXIY=-1000
      MNIX=1000
      MNIY=1000
      KNBP=0
      KNVC=0
      KNDT=0
      rktime=0.d0
C     SET TRACE-BACKUP VARIABLES
      CALL SETTRC
      MCONP=0
      POLAR=0.d0
      LPOLAR=1.d0
      CALL RMOVE(CREF,CUEF,3)
      CALL IMOVE(NRYIND,NUYIND,5)
      CALL SET2(CNVV,1.d0,20)
      DO 10 I=1,5
        PCV(I)=PFV(I)
        DO 10 J=1,4
10       CNVC(I,J)=CQVC(I,J)
C     DEFAULT WIND AND ABSORBSION IDENTIFIERS
      MODU(1)=1H  ! 10 Jun 2003
      MODU(2)=0.d0
      MODU(3)=' ' ! 2 Jun 2004
      MODU(4)=0.d0
      MODL(1)=' ' ! 2 Jun 2004
      MODL(2)=0.d0
      MODL(3)=' ' ! 2 Jun 2004
      MODL(4)=0.d0
C
C     INITIALIZE MODEL DESCRIPTIONS TO BLANKS
      CALL ISET2(IDSDC,' ' ,10) ! 10 Jun 2003
      CALL ISET2(IDSDT,' ' ,10)
      CALL ISET2(IDSDU,' ' ,10)
      CALL ISET2(IDSC,' ' ,10)
      CALL ISET2(IDSM,' ' ,10)
      CALL ISET2(IDST,' ' ,10)
      CALL ISET2(IDSU,' ' ,10)
      CALL ISET2(IDSG,' ' ,10)
      CALL ISET2(IDSDG,' ' ,10)
      CALL ISET2(IDSV,' ' ,10)
      CALL ISET2(IDSDV,' ' ,10)
      CALL ISET2(IDSPR,' ' ,10)
      CALL ISET2(IDSDP,' ' ,10)
      CALL ISET2(IDSR,' ' ,10)
      CALL ISET2(IDSS,' ' ,10)
      CALL ISET2(IDSDS,' ' ,10)
C
C     GET MACHINE CONSTANTS
C%      CALL DFCNST(RMACH,5)
C 27 Aug 2008      call mach
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY STDINI_test                                            ! 11 Mar 2009
C ----------------------------------------------------------------------!-------
      CALL OPNREP(NDEVTMP,'TAPE4')
      CALL OPNURP(NDEVBIN,'TAPE6')
      CALL OPNREP(9,'PUNCH')
C
C ----------------------------------------------------------------------!-------
      ENTRY STDINT_test                                            ! 11 Mar 2009
C ----------------------------------------------------------------------!-------
      write (*,888)
888   format (' NAME OF INPUT DATA FILE [DINP]: ',$)
      READ (*,'(A12)')dinpname
      if (dinpname.eq.'            ')dinpname='dinp'
      OPEN(UNIT=NRYIND,FILE=dinpname,STATUS='OLD',ERR=1000)
      REWIND NRYIND
C%      CALL OPNREP(2,'OUTPUT')
      CALL OPNREP(3,'DOUTP')
      CALL opnrep(NDEVGRP,'TAPE5')
C
C  INITIALIZE RAYSET FILE
      READ(NRYIND,'(A)',END=1000) INITID
C%      CALL SETW
      xyz=setw(xyz)
      zzz= SETGET(x)
      CALL SYSDAT(DAT)
      CALL SYSTIM(TOD)
      CALL SYSSEC(SEC)
C     INITIALIZE POINT MODEL LIST
      LIST(1)=1
      LIST(2)=0
C     FILL FORMAT CONTROL ARRAYS
      PNTBL(1)=1
C     ALLOW FOR AN 80 CHARACTER IDENT STRING(A8)
      PNTBL(2)=NCHPG1
C     ALLOW MAXIMUM 500 WORD PER VARIABLE
      PITBL(2)=NWPV
C     FOR BOTH X AND Y PLUS 1 FOR NUMBER OF VARIABLES
      PNTBL(3)=NSPGP
      RETURN
1000  STOP 'DINP FORMAT ERROR'
C ----------------------------------------------------------------------!-------
      END SUBROUTINE CONBLK_test                                   ! 11 Mar 2009
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE PRINTR(EVENT,CARD)
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
C     PRINTS OUTPUT AND OUTPUTS RAYSETS(MACHINE READABLE OUTPUT)
C     WHEN 'CARD' ARGUMENT NONZERO.
C ----------------------------------------------------------------------!-------
 
      implicit double precision (a-h,o-z)
 
      CHARACTER EVENT*9,NWHY*8,CC*1,PC*1,TMP*9,id*80
      LOGICAL PCNTRL                                              ! 28 Jul 2010

      DOUBLE PRECISION KNORM
      DIMENSION G0(3,3),G1(3,3)
      CHARACTER*12 HEADRS(20),HEAD(20),UNITS(20),UNIT(20),numstg
      DIMENSION RPRINT(20),NPR(20)

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/ident/id
      COMMON /RK/ NEQS,STEP,MODE,E1MAX,E1MIN,E2MAX,E2MIN,FACT,RSTART

C     COMMON DECK "CERR" INSERTED HERE
      COMMON/ERR/NERG,NERR,NERT,NERP

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 "PP" INSERTED HERE                               ! 17 Feb 2006
      DOUBLE PRECISION MODP                                        ! 17 Feb 2006
      COMMON/PP/MODP(4),Pressure,PPT,PPR,PPTH,PPPH                 ! 17 Feb 2006

C     COMMON DECK "DS" INSERTED HERE                               ! 27 Mar 2007
      DOUBLE PRECISION ds,modds                                    ! 27 Mar 2007
      COMMON/DS/MODDS(4),ds,pdst,pdsr,pdsth,pdsph                  ! 27 Mar 2007

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

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 "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

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 MMODEL,MFORM,MID
C
C     WIND        100-124
      EQUIVALENCE (W(100),UMODEL),(W(101),UFORM),(W(102),UID)
C
C     DELTA WIND        125-149
      EQUIVALENCE (W(125),DUMODEL),(W(126),DUFORM),(W(127),DUID)
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
      EQUIVALENCE(T,TPULSE),(PHI,PH),(TH,THETA)
c
      DOUBLE PRECISION gtoo(12)
      equivalence (gtoo(1),g)
 
      DATA
     2  HEADRS(7)/'  PHASE TIME'/,UNITS(7)/'   SEC'/
     3 ,HEADRS(8)/'  ABSORPTION'/,UNITS(8)/'   DB '/
     4 ,HEADRS(9)/'   DOPPLER  '/,UNITS(9)/'  C/,S '/
      DATA  HEADRS(10)/' PATH LENGTH'/,UNITS(10)/'   KM '/
     6 ,HEADRS(11)/' TERRAIN'/,UNITS(11)/' '/
     7 ,HEADRS(12)/' TERRAIN PGR'/,UNITS(12)/' '/
     8 ,HEADRS(13)/' TERRAIN PGT'/,UNITS(13)/' '/
     9 ,HEADRS(14)/' TERRAIN PGP'/,UNITS(14)/' '/
     1 ,HEADRS(15)/' LN PRESSURE'/,UNITS(15)/' '/               ! 16 Feb 2006
C
C      ROUND-OFF FUNCTION
      ROUND(X)=dSIGN(dABS(X)+0.5d0,X)
C
      NWHY=EVENT(2:)
      if (event.eq.'RRCVR    ') write (nwhy,'(a4,i4)')event(2:5),ihop
      GO TO 10
C
C     INITIALIZATION ENTRY POINT FOR PRINTR. CALLED FOR EACH NEW
C     W-ARRAY.
C ----------------------------------------------------------------------!-------
      ENTRY IPRINTR
C ----------------------------------------------------------------------!-------
C********* NEW W ARRAY -- REINITIALIZE
      PCNTRL=W(73).NE.0.d0                                        ! 28 Jul 2010
      NEWWP=.FALSE.
      SPL=dSIN (PLON-TLON)
      CPL=dSIN (PID2-(PLON-TLON))
      SP=dSIN (PLAT)
      CP=dSIN (PID2-PLAT)
      SL=dSIN (TLAT)
      CL=dSIN (PID2-TLAT)
C********* MATRIX TO ROTATE COORDINATES
      G0(1,1)=CPL*SP*CL-CP*SL
      G0(1,2)=SPL*SP
      G0(1,3)=-SL*SP*CPL-CL*CP
      G0(2,1)=-SPL*CL
      G0(2,2)=CPL
      G0(2,3)=SL*SPL
      G0(3,1)=CL*CP*CPL+SP*SL
      G0(3,2)=CP*SPL
      G0(3,3)=-SL*CP*CPL+SP*CL
      DENM=G0(1,1)*G0(2,2)*G0(3,3)+G0(1,2)*G0(3,1)*G0(2,3)
     1   +G0(2,1)*G0(3,2)*G0(1,3)-G0(2,2)*G0(3,1)*G0(1,3)
     2   -G0(1,2)*G0(2,1)*G0(3,3)-G0(1,1)*G0(3,2)*G0(2,3)
      if (denm .eq. 0.d0) stop 'denm=0. in PRINTR'
C********* THE MATRIX G1 IS THE INVERSE OF THE MATRIX G
      G1(1,1)=(G0(2,2)*G0(3,3)-G0(3,2)*G0(2,3))/DENM
      G1(1,2)=(G0(3,2)*G0(1,3)-G0(1,2)*G0(3,3))/DENM
      G1(1,3)=(G0(1,2)*G0(2,3)-G0(2,2)*G0(1,3))/DENM
      G1(2,1)=(G0(3,1)*G0(2,3)-G0(2,1)*G0(3,3))/DENM
      G1(2,2)=(G0(1,1)*G0(3,3)-G0(3,1)*G0(1,3))/DENM
      G1(2,3)=(G0(2,1)*G0(1,3)-G0(1,1)*G0(2,3))/DENM
      G1(3,1)=(G0(2,1)*G0(3,2)-G0(3,1)*G0(2,2))/DENM
      G1(3,2)=(G0(3,1)*G0(1,2)-G0(1,1)*G0(3,2))/DENM
      G1(3,3)=(G0(1,1)*G0(2,2)-G0(2,1)*G0(1,2))/DENM
      R0=EARTHR+XMTRH
C********* CARTESIAN COORDINATES OF TRANSMITTER
      XR=R0*G0(1,1)
      YR=R0*G0(2,1)
      ZR=R0*G0(3,1)
      CTHR=G0(3,1)
      STHR=dSIN (dACOS (CTHR))
      PHIR=dATAN2(YR,XR)
      ALPH=dATAN2(G0(3,2),G0(3,3))
C*********
      NR=0
      NP=0
      NERG=0
      NERR=0
      NERT=0
      NERP=0
C     INSURE NO GARBLE IN HEADERS
      HEAD(1)=' '
      UNIT(1)=' '
C
      DO 7 NN=7,20
      IF (W(NN+50).EQ.0.d0) GO TO 7
C********* DEPENDENT VARIABLE NUMBER NN IS BEING INTEGRATED
C********* NR IS THE NUMBER OF DEPENDENT VARIABLES BEING INTEGRATED
      NR=NR+1
C     ENABLE SELECTED RELATIVE ERROR PRINTOUTS FOR TERRAIN(G) OR
C     ITS DERIVATIVES WITH RESPECT TO (R)ANGE, (T)HETA OR (P)HI.
      IF(NN.EQ.11) NERG=NR
      IF(NN.EQ.12) NERR=NR
      IF(NN.EQ.13) NERT=NR
      IF(NN.EQ.14) NERP=NR
C      IF (W(NN+50).NE.2.d0) GO TO 7
      IF (W(NN+50).EQ.1.d0) GO TO 7                             ! 21 Dec 2005
C********* DEPENDENT VARIABLE NUMBER NN IS BEING INTEGRATED AND PRINTED.
C********* NP IS THE NUMBER OF DEPENDENT VARIABLES BEING INTEGRATED AND
C********* PRINTED
      NP=NP+1
C********* SAVE THE INDEX OF THE DEPENDENT VARIABLE TO PRINT
      NPR(NP)=NR
      HEAD(NP)=HEADRS(NN)
      HEAD(NP)=HEADRS(NN) ! ???
      UNIT(NP)=UNITS(NN)
    7 CONTINUE
      NPM=MIN0(NP,3)
      NP1=NPM+1
      P=0.d0
      ABSORB=0.d0
      DOPP=0.d0
      NEQS=NR+6
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY PRNHD1(CC)
C ----------------------------------------------------------------------!-------
C     PRINT PRINTR HEADER 1
C   IF NO OUTPUT NEEDED EXIT NOW
c      PRINT *, ' prtsrp = ', prtsrp
      IF(PRTSRP.NE.0.d0) RETURN
C
C     ADD NUMBER OF LINES NEEDED FOR THIS HEADER
      PC=CC
      IF(CC.EQ.'1') CALL NEWPAG(NPAG,INT(PAGLN),PC)
      LINES=LINES+3
C********* PRINT COLUMN HEADINGS
C
      if (npm.gt.0) WRITE(3,1100) PC,(HEAD(NN),NN=1,NPM)        ! 30 Apr 2006
      if (npm.le.0) WRITE(3,1100)                               ! 30 Apr 2006
C
1100  FORMAT (   A,T25,'ELEVATION',T54,'AZIMUTH',T71,'ELEVATION'/
     3 ,T20,2('ABOVE',6X),T53,'DEVIATION',T72,'ANGLE'/'  ERROR   EVENT'
     4  ,T20,'SEA LEVEL  TERRAIN     RANGE',4X,2('XMTR   LOCAL',5X),
     5  'PULSE TIME',3A12)
C
      if (npm.gt.0) WRITE(3,1150) (UNIT(NN),NN=1,NPM)           ! 30 Apr 2006
      if (npm.le.0) WRITE(3,1150) ! 30 Apr 2006
 1150 FORMAT(13X,3(8X,'KM'),2X,2(6X,'DEG',5X,'DEG'),T88
     1 ,'SEC',3X,3(3X,A7,2X))
C
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY PRNHD2(CC)
C ----------------------------------------------------------------------!-------
C     PRINT PRINTR HEADER 1
C   IF NO OUTPUT NEEDED EXIT NOW
      IF(PRTSRP.NE.0.0) RETURN
C     PAGE BY HALF LENGTH
      LINSPP=PAGLN/2
      IF(LINSPP.LT.40) LINSPP=PAGLN
      PC=CC
      IF(CC.EQ.'1') CALL NEWPAG(NPAG,LINSPP,PC)
C     ADD NUMBER OF LINES NEEDED FOR THIS SUBHEADER
      LINES=LINES+1
C
      IF(ELEND.GE.ELBEG+ELSTEP) THEN
         WRITE(3,'(A,''ELEVATION ANGLE OF TRANSMISSION =''
     1    ,F10.4,'' DEG'')') PC,BETA*DEGS  !Rev 7/16/90  remove SNGL
      ELSEIF(AZEND.GE.AZBEG+AZSTEP) THEN
         WRITE(3,'(A,''AZIMUTH ANGLE OF TRANSMISSION =''
     1    ,F10.4,'' DEG'')') PC,AZ1*DEGS   !Rev 7/16/90  remove SNGL
      ENDIF
      RETURN

C ----------------------------------------------------------------------!-------

C     IF PRINTING SUPPRESSED AND RAYSETS OFF NOTHING TO DO
10    IF(PRTSRP.NE.0.d0 .AND. CARD.EQ.0.d0) RETURN
      CALL DISPER
      IF (CARD.EQ.0.d0 .OR. IHOP.NE.0) GO TO 12
C********* OUTPUT A TRANSMITTER RAYSET
C        NOTE: THIS IS A SPECIAL CASE, ALL OTHER RAY EVENTS ARE
C           OUTPUT AT CODE 'PUNCH A RAYSET' BELOW.
C
      TLOND=TLON*DEGS
      IF (TLOND.LT.0.d0) TLOND=TLOND+360.d0
      TLATD=TLAT*DEGS
      IF (TLATD.LT.0.d0) TLATD=TLATD+360.d0
      AZ=AZ1*DEGS
      EL=BETA*DEGS
      NHOP=HOP
      NXMTRH=ROUND(XMTRH*1.d4)
      NTLATD=ROUND(TLATD*1.d3)
       NTLOND=ROUND(TLOND*1.d3)
      NRCVRH=ROUND(RCVRH*1.d4)
       NF=ROUND(OW*1.d4)
      NAZ=ROUND(AZ*1.d5)
       NEL=ROUND(EL*1.d5)
      NPOLAR1=ROUND(DREAL(POLAR)*1.D2)
       NPOLAR2=ROUND(DIMAG(POLAR)*1.D2)
      ndens = 0                                                    ! 20 Apr 2007
      if(ds.ne.0.d0) ndens = round(dlog10(ds/1.d6)*1.d5)           ! 20 Apr 2007
      WRITE(9,1201)ID(:3),TYPE(NTYP),NXMTRH,NTLATD,NTLOND,NRCVRH,NF,NAZ
     1 ,NEL,NPOLAR1,NPOLAR2,NHOP,'T'
     2 , ndens                                                     ! 20 Apr 2007
c 1201 FORMAT(A3,A1,I9,2I6,2I9,2I10,4X,I5,I4,i3,A1)
 1201 FORMAT(A3,A1,I9,2I6,2I9,2I10,4X,I5,I4,i3,A1,i10)             ! 20 Apr 2007
C*********
12    continue
!      V=1.d10
C     OBTAIN THE WORST ERROR OF THOSE ENABLED.
!      IF (KAY2.NE.0.d0) V=(KR**2+KTH**2+KPH**2)/KAY2-1.d0
        v = errRel                                                ! 27 Jul 2010
      IF(EVENT.EQ.' X wave  ') v = errRell(nsmall)  ! 21 Feb 2011 !  9 Mar 2012
      IF(EVENT.EQ.' O wave  ') v = errRell(n3)      ! 21 Feb 2011 !  9 Mar 2012

      ERT=1HK
      V=RERR(V,ERT,'G',NERG,G)
      V=RERR(V,ERT,'R',NERR,PGR)
      V=RERR(V,ERT,'T',NERT,PGTH)
      V=RERR(V,ERT,'P',NERP,PGPH)
C
      H=R-EARTHR
      STH=dSIN (THETA)
      CTH=dSIN (PID2-THETA)
C********* CARTESIAN COORDINATES OF RAY POINT, ORIGIN AT TRANSMITTER
      XP=R*STH*dSIN (PID2-PHI)-XR
      YP=R*STH*dSIN (PHI)-YR
      ZP=R*CTH-ZR
C********* CARTESIAN COORDINATES OF RAY POINT, ORIGIN AT TRANSMITTER AND
C********* ROTATED
      EPS=XP*G1(1,1)+YP*G1(1,2)+ZP*G1(1,3)
      ETA=XP*G1(2,1)+YP*G1(2,2)+ZP*G1(2,3)
      ZETA=XP*G1(3,1)+YP*G1(3,2)+ZP*G1(3,3)
      RCE2=ETA**2+ZETA**2
      RCE=dSQRT (RCE2)
C********* GROUND RANGE
      RANGE=EARTHR*dATAN2(RCE,EARTHR+EPS+XMTRH)
C********* ANGLE OF WAVE NORMAL WITH LOCAL HORIZONTAL
      ELL=dATAN2(KR,dSQRT (KTH*kth+KPH*kph))*DEGS
C********* STRAIGHT LINE DISTANCE FROM TRANSMITTER TO RAY POINT
      SR=dSQRT (RCE2+EPS*eps)
C********* TERRAIN RELATIVE HEIGHT
      GRH=GET(gtoo)/PGR
C     REPORT GROUP TIME AS FIRST 'OPTION'
      RPRINT(1)=T
! 23 Aug 2011      IF (NP1.LT.2) GO TO 16
      IF (NP.LT.1) GO TO 16                                        ! 23 Aug 2011
C     ADD MORE OPTIONS IF REQUESTED
! 23 Aug 2011      DO 15 I=2,NP1
      DO 15 I=2,NP+1                                               ! 23 Aug 2011
      NN=NPR(I-1)
   15 RPRINT(I)=RKVARS(NN)
C
c%16    IF(V.GE.0.d0) THEN
c%         WRITE(TMP,'(A1,0PE7.1)') ERT,V
c%      ELSE
16         WRITE(TMP,'(A1,0PE7.1)') ERT,V
c%      ENDIF
C     DETERMINE WHERE TO PUT A SPACE
C     AT BEGINNING OR END OF 1ST 2 VALUES
      KT=1
      IF(TMP(1:1).EQ.'K') KT=2
C
      if(prtsrp.eq.0..AND.W(73).LE.1.D0)! 22 Nov 2010 29 Nov 2010  ! 22 Feb 2011
     & write (*,222)tmp(kt:9),nwhy,h,range,beta*degs
222   format (2x,2a8,f10.4,f11.4,27x,'ELEV.= ',f8.3)               ! 23 Aug 2011

      if(PCNTRL) ks=kr*kr+kth*kth+kph*kph           ! 20 Jul 2010  ! 22 Nov 2010

      if(W(73).eq.2.D0) write (*,224)tmp(kt:9),nwhy,h,range                            
     & ,sol(n3),sol(nsmall),beta*degs               ! 22 Nov 2010  ! 29 Nov 2010

224   format (2x,2a8,f10.4,f11.4,2f13.9,1x,'ELEV.= ',f8.3)! 22 Nov 2010! 23 Aug 2011

C ----------------------------------------------------------------------!-------
      if((EVENT.EQ.' O wave  '.OR.EVENT.EQ.' X wave  ' !22 Feb 2011!  9 Mar 2012
     & .OR.EVENT.EQ.' diagnose')                                   ! 14 Mar 2011
     & .and.W(73).eq.3.D0              ! 28 Jul 2010 ! 22 Nov 2010 ! 29 Nov 2010
     & .or.W(73).eq.4.D0)                                          ! 15 Mar 2011
     & write (*,223)tmp(kt:9),nwhy,h,range,n3,n4,nbig,nsmall       ! 29 Jul 2010
     & ,Hrel
     & ,sol(1),soli(1) ,sol(2),soli(2) 
     & ,sol(3),soli(3) ,sol(4),soli(4)                             !  4 Aug 2010
     & , dd(0)
     & ,solb(1),solbi(1) ,solb(2),solbi(2) 
     & ,solb(3),solbi(3) ,solb(4),solbi(4)                         !  4 Aug 2010
     & , dd(1)
     & ,sola(1),solai(1) ,sola(2),solai(2) 
     & ,sola(3),solai(3) ,sola(4),solai(4)                         !  4 Aug 2010
     & ,dd(2) 
     & ,kaysq(1),kaysqi(1) ,kaysq(2),kaysqi(2) 
     & ,kaysq(3),kaysqi(3) ,kaysq(4),kaysqi(4) 
     & ,dd(3) 
     & ,kaysqb(1),kaysqbi(1) ,kaysqb(2),kaysqbi(2) 
     & ,kaysqb(3),kaysqbi(3) ,kaysqb(4),kaysqbi(4) 
     & ,dd(4) 
     & ,kaysqa(1),kaysqai(1) ,kaysqa(2),kaysqai(2) 
     & ,kaysqa(3),kaysqai(3) ,kaysqa(4),kaysqai(4) 
     & ,alpha
     & ,relerr 
     & ,ray 
     & ,relerrb 
     & ,ray 
     & ,relerra                                                    ! 29 Jul 2010
     & ,ray 
     & ,errRell
     & ,ray 
     & ,errRelb 
     & ,ray 
     & ,errRela                                                    ! 19 Nov 2010
     & ,ray 
     & ,ks 
     & ,kay2                                                       ! 29 Jul 2010
!      if(W(73).eq.4.d0) print*,'root',int(ray),
!     & 'KR = ',KR, 'KTH = ',KTH, 'KPH = ',KPH                      ! 17 Jul 2018
!     & ,'r=',r,'th=',th,'ph=',ph,'kay2=',kay2                      ! 17 Jul 2018
223   format (1x,2a8,f10.4,f11.4
     & ,', n3=',i2,', n4=',i2,', nbig=',i2,', nsmall=',i2 
     & /(1x,0PE7.1,8g9.2))                                         ! 29 Jul 2010
      if(EVENT.EQ.' O wave  '.and.prtsrp.eq.0..and.W(73).eq.3.d0)  !  9 Mar 2012
     & write (*,223)                                ! 22 Nov 2010  ! 29 Nov 2010

C     TEST IF NEED NEW PAGE
      CALL TSTPAG(3)
      IF (SR.GE.1.d-6) GO TO 20
C********* TOO CLOSE TO TRANSMITTER TO CALCULATE DIRECTION FROM
C********* TRANSMITTER
      IF(PRTSRP.EQ.0.d0)        !Rev 7/16/90 remove SNGL below
     1 WRITE(3,1500)TMP(KT:9),NWHY,H,GRH,RANGE,ELL,(RPRINT(NN),NN=1,NP1)
 1500 FORMAT (1X,A8,A8,2F10.4,F11.4,26X,F8.3,4F12.4)
C     SET RAYSET VARIABLES TO UNDEFINDED VALUES FOR FLAGS
      AZDEV=999.d0
      AZA=999.d0
      GO TO 40
C********* ELEVATION ANGLE OF RAY POINT FROM TRANSMITTER
   20 EL=dATAN2(EPS,RCE)*DEGS
      IF (RCE.GE.1.d-6) GO TO 30
C********* NEARLY DIRECTLY ABOVE OR BELOW TRANSMITTER. CAN NOT CALCULATE
C********* AZIMUTH DIRECTION FROM TRANSMITTER ACCURATELY
      IF(PRTSRP.EQ.0.d0)    !Rev 7/16/90 remove SNGL below      
     1   WRITE(3,2500)TMP(KT:9),NWHY,H,GRH,RANGE,EL,ELL,(RPRINT(NN),
     2   NN=1,NP1)
 2500 FORMAT (1X,A8,A8,2F10.4,F11.4,17X,F9.3,F8.3,
     1 4F12.4)
      GO TO 40
C********* AZIMUTH ANGLE OF RAY POINT FROM TRANSMITTER
   30 ANGA=dATAN2(ETA,ZETA)
      AZDEV=180.d0-DMOD(540.d0-(AZ1-ANGA)*DEGS,360.D0)
      IF (KTH.NE.0.d0.OR.KPH.NE.0.d0) GO TO 34
C********* WAVE NORMAL IS VERTICAL, SO AZIMUTH DIRECTION CANNOT BE
C********* CALCULATED
      IF(PRTSRP.EQ.0.d0)    !Rev 7/16/90 remove SNGL below
     1 WRITE(3,3000)TMP(KT:9),NWHY,H,GRH,RANGE,
     1 AZDEV,EL,ELL,(RPRINT(NN),NN=1,NP1)
 3000 FORMAT (1X,A8,A8,2F10.4,F11.4,F9.3,8X,F8.3,F8.3,
     1 4F12.4)
      GO TO 40
   34 ANA=ANGA-ALPH
      SANA=dSIN (ANA)
      SPHI=SANA*STHR/STH
      CPHI=-dSIN (PID2-ANA)*dSIN (PID2-(PHI-PHIR))+SANA*dSIN (PHI-PHIR)
     1 *CTHR
      AZA=180.d0-DMOD(540.d0-(DATAN2(SPHI,CPHI)-DATAN2(KPH,KTH))*DEGS,
     1 360.D0)
      IF(PRTSRP.EQ.0.d0)   !Rev 7/16/90 remove SNGL below
     1  WRITE(3,3500)TMP(KT:9),NWHY,H,GRH,RANGE,
     1  AZDEV,AZA,EL,ELL
     2  ,(RPRINT(NN),NN=1,NP1)
 3500 FORMAT (1X,A8,A8,2F10.4,F11.4,2(F9.3,F8.3),
     1 4F12.4)
C*********
   40 LINES=LINES+1
      IF (NP.LE.3) GO TO 45
C********* ADDITIONAL LINE TO PRINT REMAINING DEPENDENT INTEGRATION
C********* VARIABLES
      IF(PRTSRP.EQ.0.d0)   !Rev 7/16/90 remove SNGL below
     1   WRITE(3,4000)(RPRINT(NN),NN=5,NP+1)                       ! 23 Aug 2011
! 23 Aug 2011     1   WRITE(3,4000)(RPRINT(NN),NN=4,NP)
 4000 FORMAT (99X,3F12.4)
      LINES=LINES+1
C
C     IF NO 'CARDS' WANTED OR AT TRANSMITTER, NO RAYSET OUTPUT
C
45    IF (CARD.EQ.0.d0 .OR. IHOP.LT.1) RETURN
C

      IF (W(65).NE.0.D0) write(3,1700) pressure                    ! 21 Feb 2006
 1700 FORMAT('Pressure = ',e25.10,' Pa')                           ! 17 Feb 2006

      IF (W(65).NE.0.D0) write(3,1800) ds/1.D6                     ! 27 Mar 2007
 1800 FORMAT('Density = ',e25.10,' kg/m^3')                        ! 27 Mar 2007

      IF (W(65).NE.0.D0) write(3,1900) dsqrt(cs)                   ! 15 May 2007
 1900 FORMAT('Sound speed = ',e25.10,' km/s')                      ! 15 May 2007

C********* PUNCH A RAYSET
      IF (AZDEV.LT.-90.d0) AZDEV=AZDEV+360.d0
      IF (AZA.LT.-90.d0) AZA=AZA+360.d0
      NR=0
      IF (W(57).EQ.0.d0) GO TO 47
C********* PHASE PATH
      NR=NR+1
      P=RKVARS(NR)
   47 IF (W(58).EQ.0.d0) GO TO 48
C********* ABSORPTION
      NR=NR+1
      ABSORB=RKVARS(NR)
C********* DOPPLER SHIFT
   48 IF (W(59).NE.0.d0) DOPP=RKVARS(NR+1)
      NHPUNCH=ROUND(HPUNCH*1.d4)
       NRANGE=ROUND(RANGE*1.d4)
      NAZDEV=ROUND(AZDEV*1.d3)
       NAZA=ROUND(AZA*1.d3)
       NELL=ROUND(ELL*1.d3)
      IF(NWHY.EQ.'GRND REF') NELL =
     1 ROUND((PID2 -DACOS(KNORM(G)/DSQRT(KR*KR+KTH*KTH+KPH*KPH)))
     2   *DEGS*1.d3)
      NABSORB=DMIN1(999999.D0,ROUND(ABSORB*1.D3))
      NDOPP=ROUND(DOPP*1.D3)
      NPOLAR1=ROUND(DREAL(POLAR)*1.D2)
      NPOLAR2=ROUND(DIMAG(POLAR)*1.D2)
      JP=ROUND(P*1.D5)
      JT=ROUND(T*1.D5)
      ndens = 0                                                    ! 20 Apr 2007
      if(ds.ne.0.d0) ndens = round(dlog10(ds/1.d6)*1.d5)           ! 20 Apr 2007
      WRITE(9,4501) NHPUNCH,NRANGE,NAZDEV,NAZA,NELL,JT,JP,NABSORB
     1 ,NDOPP,NPOLAR1,NPOLAR2,IHOP,EVENT(1:1)
     2 , ndens                                                     ! 20 Apr 2007
c 4501 FORMAT(2I9,3I6,2I10,2I6,2i4,i3,A1)
 4501 FORMAT(2I9,3I6,2I10,2I6,2i4,i3,A1,i10)                       ! 20 Apr 2007
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PRINTR
C ----------------------------------------------------------------------!-------
 
 
