C file S0.FOR
C**********************************************************************
! 29 Jun 2010 - Go to next elevation angle in evanescent region
C 30 Apr 2009 - replaced REAL*8 with DOUBLE PRECISION
C 24 Mar 2008 - Initialize routine SREAD1
C 13 Mar 2008 - Corrected length in common block /PCONST/
C  7 Dec 2006 - Added initializing of path length and terrain function
C  7 Dec 2006 - Corrected initializing of dependent variables
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 DFSYS, READW1, GTUNIT, SREAD1, READW to S5.for
C**************************************************************************** 
Contains  RAYTRC
C**********************************************************************
C 17 Feb 2006 - Added w(556)=nv
C 16 Feb 2006 - Added initialization for integration of pressure
C 31 January 2006 - Added /B13/, /B14/, /B15/, and /B16/ to READW
C 5 January 2006 - Changed printing of CONV to double precision in READW
C 21 December 2005 - Initializing variables to coincide with EIGEN
C  6 Sept 2005 - The product of W42 (max single step error) and
C                W46 (min integration step size) has to be greater than
C                10E4 (Earth radius) / 10E15 (15-digit double precision)
C                This equals 10E-11 for double precision in FORTRAN 90
C 28 Nov 2003 (marked) corrected relative-to-terrain transmitter heights
C (SUBROUTINE RAYTRC and FUNCTION READW)
C 10 Jun 2003 (not marked) changed lines 381 (SUBROUTINE DFSYS)
C 10 Jun 2003 (not marked) changed lines 927, 928, 935 (FUNCTION READW)
C**********************************************************************
      PROGRAM RAYTRC
C**********************************************************************
C     MAIN PROGRAM FOR THE RAY TRACING PACKAGE.
C          SETS THE INITIAL CONDITIONS FOR EACH RAY AND CALLS TRACE
C
C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
      IMPLICIT double precision (a-h,o-z)

      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN
      LOGICAL WCHANGE,FIRST
      REAL(KIND=8) :: WS(400),kay2,kay2i,kr,kth,kph 
     & , kz2a,kz2ai,kzb,kzbi        ! 16 Mar 2006
     1,maxstp,maxerr,intyp,llat,llon
      DOUBLE PRECISION modg
C     COMMON DECK "FILEC" INSERTED HERE
      COMMON/FILEC/NPLTDP
C     COMMON DECK "CERR" INSERTED HERE
      COMMON/ERR/NERG,NERR,NERT,NERP

C     COMMON DECK "GG" INSERTED HERE
      COMMON/GG/MODG(4)
      COMMON/GG/G,PGR,PGRR,PGRTH,PGRPH
      COMMON/GG/PGTH,PGPH,PGTHTH,PGPHPH,PGTHPH,GSELECT,GTIME

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

C     COMMON DECK "HDR" INSERTED HERE
      CHARACTER*10 INITID*80,DAT,TOD
      COMMON/HDR/SEC,stsec                             !Rev 7/13/90
      COMMON/HDRC/INITID,DAT,TOD

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

C     COMMON DECK "WWR" INSERTED HERE
      COMMON/WW/DUM(10),MAXW,W(1000)
      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))
      EQUIVALENCE (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
      CHARACTER TMP80*80
C
      DATA FIRST/.TRUE./
C
C   SET CONSTANTS , INITITALIZE FILES  AND MAKE SYSTEM CALLS
c%      CALL SYSINI
      CALL CONBLK
      CALL SETDSP                   ! Initialize dispersion relation routine
      CALL SETSRD ! Initialize routine SREAD1                      ! 24 Mar 2008
      CALL STDINI
      call syssec(stsec)               !Rev 7/13/90  get start time
C     INITIALIZE LINE COUNTS FOR CURRENT AND LAST HEADER
      LINES=0
      LHDRPG=0
C
      REWIND NRYIND
C  INITIALIZE RAYSET FILE
5     READ(NRYIND,'(A)',END=8) TMP80
      WRITE(9,'(A)') TMP80
      GO TO 5
C
8     WRITE(9,'(A)') '****SOFT EOF***+++++'
      REWIND NRYIND
C     READ IN USERS NAME AND TELEPHONE EXTENSION FOR IDENTIFYING
C     MICROFILM PLOTS.
C
C     REPOSITION DINP AFTER IDENTIFIER
      READ(NRYIND,'(A)',END=6000) TMP80
C
C********* READ W ARRAY AND PRINT NON-ZERO VALUES
      IRUN=0
      NSET=0
C
10    xyz= READW(xx)
      NSET=NSET+1
      ICODE=ND2B(IDINT(RAYFNC))
C PROCESS RAYPATH CALCULATIONS ONLY IF W(29) IS BEING USED
c*      IF((RAYFNC.NE.0.).AND.(AND(ICODE,4).EQ.0)) GO TO 10
C*  For VAX A Newhall  (AND to IAND)
      IF((RAYFNC.NE.0.d0).AND.(IAND(ICODE,4).EQ.0)) GO TO 10
C
C********* LET ROUTINES PRINTR AND RAYPLT KNOW THERE IS A NEW W ARRAY
      NEWWP=.TRUE.
      NEWWR=.TRUE.
      IRUN=IRUN+1
C
C     SET BINARY RAY FILE UNIT TO ZERO IF NO OUTPUT WANTED
      NPLTDP=0
      IF(BINRAY.NE.0.d0) NPLTDP=NDEVBIN
C
C******** INITIALIZE THE MODELS VIA 'DISPER'
      CALL IDISPER
C
      IF(FIRST.AND.PLT.NE.0.d0) CALL DDINIT(8,INITID)
      IF(PLT.NE.0.d0) FIRST=.FALSE.
C
      OW=0.d0
      BETA=0.d0
      AZ1=0.d0
C
      CALL HEADER1
C
C********* PRINT OUT THE CONTENTS OF THE 'W' ARRAY
C********* DETERMINE TRANSMITTER LOCATION IN COMPUTATIONAL COORDINATE
C********* SYSTEM (GEOMAGNETIC COORDINATES IF DIPOLE FIELD IS USED)
      SP=DSIN (PLAT)
      CP=DSIN (PID2-PLAT)
      SDPH=DSIN (TLON-PLON)
      CDPH=DSIN (PID2-(TLON-PLON))
      SL=DSIN (TLAT)
      CL=DSIN (PID2-TLAT)
      ALPHA=DATAN2(-SDPH*CP,-CDPH*CP*SL+SP*CL)
      TH0=DACOS (CDPH*CP*CL+SP*SL)
      PH0=DATAN2(SDPH*CL,CDPH*SP*CL-CP*SL)
C
      R=EARTHR
      TH=TH0
      PH=PH0
      CALL TOPOG
C
C     OBTAIN ABSOLUTE HEIGHT OF THE TRANSMITTER.
C     IF IT WAS SPECIFIED AS RELATIVE TO THE TERRAIN, THEN REMOVE
C**   THE FLAG VALUE 10**40 WHICH WAS ADDED AT INPUT.
C     THE FLAG VALUE 10**10 WHICH WAS ADDED AT INPUT.
C
      TMP=XMTRH
C**   IF(XMTRH .EQ. 1.E-40) XMTRH=0.0
c     IF(XMTRH .EQ. 1.d-10) XMTRH=0.0
      IF(XMTRH .EQ. 1.d-25) XMTRH=0.0                ! 28 Nov 2003
C**   IF(DABS(XMTRH) .GE. 1.E20 ) XMTRH=XMTRH*1.E-40
c     IF(DABS(XMTRH) .GE. 1.d10 ) XMTRH=XMTRH*1.d-10
      IF(DABS(XMTRH) .GE. 1.d10 ) XMTRH=XMTRH*1.d-25 ! 28 Nov 2003
      IF(TMP.NE.XMTRH) XMTRH=XMTRH-G/PGR
C
C     CHECK THAT TRANSMITTER IS ABOVE TERRAIN.
C
      IF(-G/PGR .LE. XMTRH) GO TO 655
C
      PRINT *,'transmitter below terrain; run #',irun
      WRITE(3,640) IRUN
640   FORMAT('0****** TRANSMITTER BELOW TERRAIN. RUN ',I3
     1 ,' TERMINATED.'/'SEE W-ARRAY PRINTOUT. INPUT CONTINUES.'//)
      xyy= setovr(xx)
      zyx= printw(xx)
      GO TO 10
C
C     PREVENT RERUNNING SAME MODEL CASES
655   IF(IRUN.EQ.1) GO TO 12
      IF(RUNSUP.NE.0.d0 .AND. .NOT.WCHANGE(WS(1),W(1))) THEN
         PRINT *,'SUPPRESSING REDUNDANT RUNSET #',NSET
         GO TO 10
      ENDIF
C
12    PRINT *,'PROCESSING FOR RUN SET NUMBER',NSET
C
      CALL RMOVE(WS,W,400)
C
      xyz=setovr (xx)
      zyx= printw(xx)
      CALL IPRINTR
C******** INITIALIZE PRINT CONTROL PARAMTERS
      LINSPP=PAGLN
      LNPHPG=LINSPP/2
      IF(LNPHPG.LT.40) LNPHPG=LINSPP
C
C********* LOOP ON FREQUENCY, AZIMUTH ANGLE, AND ELEVATION ANGLE
      NFREQ=1
      IF (FSTEP.NE.0.d0) NFREQ=(FEND-FBEG)/FSTEP+1.5d0
      NAZ=1
      IF (AZSTEP.NE.0.d0) NAZ=(AZEND-AZBEG)/AZSTEP+1.5d0
      NBETA=1
      IF (ELSTEP.NE.0.d0) NBETA=(ELEND-ELBEG)/ELSTEP+1.5d0
      DO 50 NF=1,NFREQ
      OW=FBEG+(NF-1)*FSTEP
      DO 45 J=1,NAZ
      AZ1=AZBEG+(J-1)*AZSTEP
      AZA=AZ1*DEGS
      GAMMA1=PI-AZ1+ALPHA
      SGAMMA=DSIN (GAMMA1)
      CGAMMA=DSIN (PID2-GAMMA1)
      DO 40 I=1,NBETA
      BETA=ELBEG+(I-1)*ELSTEP
      EL=BETA*DEGS
      print *,'ELEVATION ANGLE=',sngl(EL),'  DEG'  !Rev 7/16/90
      CBETA=DSIN (PID2-BETA)
      R=EARTHR+XMTRH
      TH=TH0
      PH=PH0
      KR=DSIN (BETA)
c      print *,'beta = ',beta
c      print *,'pid2 = ',pid2
c      print *,'cbeta = ',cbeta
      KTH=CBETA*CGAMMA
      KPH=CBETA*SGAMMA
      TPULSE=0.d0
      RSTART=1.d0
C********* THE FOLLOWING LINE NEEDED FOR RAY TRACING IN COMPLEX SPACE
      SGN=1.d0
C********* CALL MODELS
!      print *,'before disper '
!      print *, ' kay2 = ', kay2
      CALL DISPER
!      print *, ' kay2 = ', kay2
!      print *,'after disper '

C...............LINPG= # of lines printed so far on this page
C...............LHDRPG= # of pages times lines per page
C...............LINES= total # of lines printed, assuming full pages
C...................(except for this page)
C...............LINSPP= # of lines per page
C...............LNPHPG= # of lines per half page
      lhdrpg=lines/linspp*linspp ! added to get correct line count
      LINPG=LINES-LHDRPG
      IF(I.EQ.1  .OR.  LINPG.GE.LINSPP-20) THEN
C
C     PUT OUT SUBHEADERS FROM MEDIA AND PRINTR ROUTINES
c      print *, 'before header2'
         CALL HEADER2
c      print *, 'after header2'
c      print *, 'before prnhd1'
         CALL PRNHD1(' ')
c      print *, 'after prnhd1'
C
C     COMPUTE LINE COUNT OF THIS HEADER
         LHDRPG=LINES/LINSPP*LINSPP
      ELSEIF(LINPG.GE.LNPHPG-10 .AND. LINPG.LE.LNPHPG) THEN
C     PUT OUT PAGE FEED WITH SUBHEADER IF AT HALF PAGE
c      print *, ' before prnhd2'
         CALL PRNHD2('1')
c      print *, ' after prnhd2'
      ELSE
C     PUT OUT SUBHEADER
c      print *, ' before prnhd2 - second'
         CALL PRNHD2(' ')
c      print *, ' after prnhd2 - second'
      ENDIF
C
      print *, ' kay2 = ', kay2
      IF (KAY2.GT.0.d0) GO TO 30
      WRITE(3,2900) OMEGMIN,OMEGMAX
c     1  ,t
 2900 FORMAT (58H0TRANSMITTER IN EVANESCENT REGION, TRANSMISSION IMPOSSI
     1BLE/20H0MINIMUM FREQUENCY =,E17.10,20H MAXIMUM FREQUENCY =,E17.10
c     2 , ' temperature = ', E17.10
     3 )
! 29 Jun 2010      GO TO 44
      GO TO 40                                                    ! 29 Jun 2010
30    r999=renorm(kr,kay2,3)
c30    CALL RENORM(KR,KAY2,3)
c      PRINT *, ' neqs = ', neqs
      if (neqs-6 .gt. 0) then
        CALL CLEAR(RKVARS,NEQS-6)
!21 December 2005 - Initializing variables to coincide with EIGEN
C********** GROUP TIME
        IF (W(56).ne.0.d0) then                      ! 7 Dec 2006
           if(W(56).ne.1.d0) TPULSE=W(56)            ! 7 Dec 2006
        ENDIF
        NV=0
C********** PHASE PATH
        IF (W(57).ne.0.d0) then                      ! 7 Dec 2006
           NV=NV+1
           if(W(57).ne.1.d0) RKVARS(NV)=W(57)        ! 7 Dec 2006
        ENDIF
C********** ABSORPTION
        IF (W(58).ne.0.d0) then                      ! 7 Dec 2006
           NV=NV+1
           if(W(58).ne.1.d0) RKVARS(NV)=W(58)        ! 7 Dec 2006
        ENDIF
C********** DOPPLER SHIFT
        IF (W(59).ne.0.d0) then                      ! 7 Dec 2006
           NV=NV+1
           if(W(59).ne.1.d0) RKVARS(NV)=W(59)        ! 7 Dec 2006
        ENDIF
C********** PATH LENGTH                              ! 7 Dec 2006
        IF (W(60).ne.0.d0) then
           NV=NV+1
           if(W(60).ne.1.d0) RKVARS(NV)=W(60)
        ENDIF
C********** TERRAIN FUNCTION                         ! 7 Dec 2006
        IF (W(61).ne.0.d0) then
           NV=NV+1
           if(W(61).ne.1.d0) RKVARS(NV)=W(61)
        ENDIF
C********** TERRAIN FUNCTION                         ! 7 Dec 2006
        IF (W(62).ne.0.d0) then
           NV=NV+1
           if(W(62).ne.1.d0) RKVARS(NV)=W(62)
        ENDIF
C********** TERRAIN FUNCTION                         ! 7 Dec 2006
        IF (W(63).ne.0.d0) then
           NV=NV+1
           if(W(63).ne.1.d0) RKVARS(NV)=W(63)
        ENDIF
C********** TERRAIN FUNCTION                         ! 7 Dec 2006
        IF (W(64).ne.0.d0) then
           NV=NV+1
           if(W(64).ne.1.d0) RKVARS(NV)=W(64)
        ENDIF
C********** Pressure ! 16 Feb 2006
        IF (W(65).ne.0.d0) then                      ! 7 Dec 2006
           NV=NV+1 ! 16 Feb 2006
           w(556) = nv ! 17 Feb 2006
           if(W(65).ne.1.d0) RKVARS(NV)=W(65)        ! 7 Dec 2006
        ENDIF
      endif
      CALL TOPOG
      IF(NERG.GT.0) RKVARS(NERG)=G
      IF(NERR.GT.0) RKVARS(NERR)=PGR
      IF(NERT.GT.0) RKVARS(NERT)=PGTH
      IF(NERP.GT.0) RKVARS(NERP)=PGPH
C
C     CALCULATE ONE RAY PATH
      CALL TRACE
      OSEC=SEC
      CALL SYSSEC(SEC)
      DIFF=SEC-OSEC
C
C     ADD TO LINES COUNT FOR ELAPSED TIME REPORT
      LINES=LINES+2
      WRITE(3,3500)SNGL(DIFF)
3500  FORMAT (/T93,'THIS RAY CALCULATION TOOK ',F8.3,' SEC')
C
      IF (PENET.AND.ONLY.NE.0.d0.AND.IHOP.EQ.1) GO TO 44
   40 CONTINUE
44    IF(PLT.GT.0.d0.AND.(NAZ.LE.1.OR.NBETA.GT.1)) CALL ENDPLT
   45 CONTINUE
      IF(PLT.GT.0.d0.AND.NAZ.GT.1.AND.NBETA.LE.1) CALL ENDPLT
      IF(PENET.AND.ONLY.NE.0.d0.AND.IHOP.EQ.1.AND.NAZ.EQ.1.AND.
     1 NBETA.EQ.1) GO TO 55
   50 CONTINUE
  55  IF (RAYSET.NE.0.d0) WRITE(9,5000)
 5000 FORMAT (78X,1H-)
      GO TO 10
6000  PRINT *,'DINP EMPTY OR NOT FOUND'
      STOP
      END PROGRAM RAYTRC
 
 
 
 
 
