C file PCPROF1.FOR
C ----------------------------------------------------------------------!-------
C Includes DFCNST
C ----------------------------------------------------------------------!-------
C 10 Nov 2009 - Added:     double precision d1mach
C  4 May 2009 - replaced REAL*8 with DOUBLE PRECISION
C 26 Aug 2008 - Replaced using rmach() with direct calls to 
C                I1mach(), R1mach(), or D1mach()
C 26 Aug 2008 - Correctly calculated machine constants using 
C                I1mach.for, R1mach.for, and D1mach.for
C 13 Mar 2008 - Corrected length in common block /TRAC/
C 24 Mar 2006 - Moved DDEND and DDTEXT  to S6
C 24 Mar 2006 - Moved C to S7
C 24 Mar 2006 - Moved COMPRS, CURVE, DONEPL, ENDPL, MX1ALF, MX2ALF, NCLOSE, 
C               NOBRDR, PAGE, PHYSOR to S8
C 24 Mar 2006 - Changed COMMON block /RAYDEV/
C 20 Mar 2006 - Altered SUBROUTINE DDTEXT
C  2 Mar 2006
C  1 Mar 2006 - Removed SETTRC, GET1
C 27 Feb 2006 - Removed DFSYS, CLEAR, ND2B
C 27 Feb 2006 - Removed GAUSEL, NUMSTG, RERROR, OPNREP, OPNURP, ZAPFIL, OVERRD
C 3 Jun 2004 - Added FUNCTION RERROR from S3.FOR
C 2 Jun 2004 - Added FUNCTION GET1 from S2.FOR
C 3 May 2004 - /CRMACH/RMACH(5) --> (9)
C ----------------------------------------------------------------------!-------
      SUBROUTINE CLEAR_TEST(A,N) ! 27 Feb 2006
C ----------------------------------------------------------------------!-------
C     SET N ELEMENTS OF ARRAY A TO ZERO
      DOUBLE PRECISION A(N)
      IF(N.LE.0) RETURN
      DO    10 I=1,N
10    A(I)=0.0
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE CLEAR_TEST
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE SETTRC_test
C ----------------------------------------------------------------------!-------
C     INITIALIZE /TRAC/ LABELED COMMON
C     COMMON DECK "TRAC" INSERTED HERE
      implicit double precision (a-h,o-z)                          ! 13 Mar 2008
      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
      GROUND=.FALSE.
      SURF=.FALSE.
C ----------------------------------------------------------------------!-------
      END SUBROUTINE SETTRC_test
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE DDEND_test
C ----------------------------------------------------------------------!-------
C     EMPTIES PLOT BUFFER AND RELEASES PLOTTING COMMAND FILE(DDPLOT)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C
c     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN,NDEVTXT ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      COMMON/KNKN/KNBP,KNVC,KNDT
      COMMON/DDLIM/MXIX,MXIY,MNIX,MNIY
C
      WRITE(NDEVTXT,100) -1
100   FORMAT(I2)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE DDEND_test
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE DDTEXT_test(N,TEXT)
C ----------------------------------------------------------------------!-------
C     WRITES A CHARACTER ARRAY PACKED A10(DDPLOT)
C     WRITE PARAMETERS TO GRAPHICS CALLS FILE
C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
 
      implicit double precision (a-h,o-z)
 
c     COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN,NDEVTXT ! 24 Mar 2006
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN         ! 24 Mar 2006
      CHARACTER TEXT*(*)
      COMMON/DD/IN,IOR,IT,IS,IC,ICC,IX,IY
      COMMON/KNKN/KNBP,KNVC,KNDT
100   FORMAT(A10,2I5,I2,1X,8A10)
      KNDT=KNDT+1
c%      WRITE(ndevtxt) 3,IX,IY

c      IX = 0                                  ! 20 Mar 2006
c      IY = 0                                  ! 20 Mar 2006
c      IOR = 0                                 ! 20 Mar 2006

      write (ndevtxt,33)3,ix,iy
33    format (i2,2i4)
      LN=MIN0(N*10,LEN(TEXT))
      IF(LN.EQ.N*10) THEN
c%         WRITE(NDEVGRP) IOR,N,N*10,(TEXT(I:I),I=1,LN)
      write (ndevtxt,34)ior,n*10,(text(i:i),i=1,ln)
34    format (2i3/80a1)
      ELSE
C%         WRITE(NDEVGRP)IOR,N,N*10,(TEXT(I:I),I=1,LN)
c%     1   ,(' ',I=LN+1,N*10)
      write (ndevtxt,35)ior,n*10,(text(i:i),i=1,ln)
     1     ,(' ',i=ln+1,n*10)
35    format (2i3/80a1/80a1)
      ENDIF
C ----------------------------------------------------------------------!-------
      END SUBROUTINE DDTEXT_test
C ----------------------------------------------------------------------!-------

C ----------------------------------------------------------------------!-------
      SUBROUTINE DFCNST(DRMACH,MAX)
C ----------------------------------------------------------------------!-------
      DOUBLE PRECISION DRMACH(MAX)
C     RETURN MAX MACHINE CONSTANTS(OR THE NUMBER AVAILABLE)
C     IN ARRAY DRMACH.

C     THIS MODULE PROVIDES CONSTANTS WHICH ARE EXPLICITLY DEPENDENT
C     ON MACHINE ARCHITECTURE AND SHOULD BE EXAMINED WHEN MAKING ANY
C     MIGRATION OF THE RAY TRACING ROUTINES.
C
C
      double precision d1mach                                     ! 10 Nov 2009
C      DOUBLE PRECISION RMACH(5) ! 3 May 2004
      DOUBLE PRECISION RMACH(9) ! 3 May 2004
      common /crmach/rmach
C
C     SINGLE PRECISION MACHINE CONSTANTS
C      (SEE DIGITAL SIGNAL PROCESSING, IEEE PRESS 1979 P S-7)
C
C     THE CONSTANTS ARE DERIVED FROM THE FOLLOWING
C     FORM FOR FLOATING POINT NUMBERS: 
C
C        SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C     FOR MOST MACHINES THE BASE B=2.
C     FOR THE CYBER 60-BIT MACHINES
C        T=48, EMIN=-974 AND EMAX=1070.
C     FOR THE PDP-11 AND VAX-11 MACHINES
C        T=24, EMIN=-127 AND EMAX=127
C
C     RMACH(1)  SMALLEST POSITIVE MAGNITUDE = B**(EMIN-1)
C     RMACH(2)  LARGEST MAGNITUDE           = B**EMAX*(1-B**(-T))
C     RMACH(3)  SMALLEST RELATIVE SPACING   = B**(-T)
C     RMACH(4)  LARGEST RELATIVE SPACING    = B**(1-T)
C     RMACH(5)                              = LOG10(B=2)
C
C     CONSTANTS HERE ARE FOR CDC 6000/7000/CYBER SERIES COMPUTERS
C     PLEASE SUBSTITUTE FOR OTHER MACHINES
C
C
C      INTEGER SMALL(2)
C      INTEGER LARGE(2)
C      INTEGER RIGHT(2)
C      INTEGER DIVER(2)
C      INTEGER LOG10(2)
CC
C      EQUIVALENCE (RMACH(1),SMALL(1))
C      EQUIVALENCE (RMACH(2),LARGE(1))
C      EQUIVALENCE (RMACH(3),RIGHT(1))
C      EQUIVALENCE (RMACH(4),DIVER(1))
C      EQUIVALENCE (RMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR CYBER 170/180 COMPUTERS
C
C*****      DATA RMACH(1) / O"00014000000000000000" /
C*****      DATA RMACH(2) / O"37767777777777777777" /
C*****      DATA RMACH(3) / O"16404000000000000000" /
C*****      DATA RMACH(4) / O"16414000000000000000" /
C*****      DATA RMACH(5) / O"17164642023241175720" /
C
C***REFERENCES  FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR
C                 A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-
C                 MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,
C                 PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  R1MACH
C     MACHINE CONSTANTS FOR VAX 11/780
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1), SMALL(2) /        128,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /
C     DATA DIVER(1), DIVER(2) /       9472,           0 /
C     DATA LOG10(1), LOG10(2) /  546979738,  -805796613 /
C
C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB /
C
C   MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING)
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1), SMALL(2) /         16,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /      15552,           0 /
C     DATA DIVER(1), DIVER(2) /      15568,           0 /
C     DATA LOG10(1), LOG10(2) /  1142112243, 2046775455 /
C
C     DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F /
C
C     MACHINE CONSTANTS FOR THE CRAY 1
C
C     DATA SMALL(1) / 201354000000000000000B /
C     DATA SMALL(2) / 000000000000000000000B /
C
C     DATA LARGE(1) / 577767777777777777777B /
C     DATA LARGE(2) / 000007777777777777774B /
C
C     DATA RIGHT(1) / 376434000000000000000B /
C     DATA RIGHT(2) / 000000000000000000000B /
C
C     DATA DIVER(1) / 376444000000000000000B /
C     DATA DIVER(2) / 000000000000000000000B /
C
C     DATA LOG10(1) / 377774642023241175717B /
C     DATA LOG10(2) / 000007571421742254654B /
C
C     The only time this is used is by GET function which is not even
C     called by the PROFIL program.  This RMACH value is included
C     just for "completeness".

c.....smallest relative spacing
C 26 Aug 2008      rmach(3) = 2.d0**(-52)
      rmach(3) = d1mach(3)                                         ! 26 Aug 2008

C     The following is equivalenced to BIGVAL in READW function, although
C     it is also probably never going to be used by the PROFIL program.

C.....largest 4-byte integer
C 26 Aug 2008      rmach(5) = 2147483647.d0
      rmach(6) = I1mach(9)                                         ! 26 Aug 2008

      CALL RMOVE(DRMACH,RMACH,MIN0(MAX,5))
C
C ----------------------------------------------------------------------!-------
      END SUBROUTINE DFCNST
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE DFSYS_test ! 27 fEB 2006
C ----------------------------------------------------------------------!-------
C*****CHANGED THE CHARACTER DECLARATION SO IT WOULD COMPILE. JAH.
C*****
C     This entire routine has been commented out except for the SYSINI
C     entry point (which actually does nothing but is there in case
C     any special system calls are required.	Most of these
C     functions have no importance for PROFIL calculations or plotting.
C
c*****	CHARACTER*(*) DAT, TIM*(*)
c*****	CHARACTER DATX*10
c*****	CHARACTER TIMX*10
C
C     ALL ROUTINES WHICH RELY ON SYSTEM BASED FUNCTIONS, SUCH AS TIME
C     OF DAY OR MAY BENEFIT FROM SYSTEM AVAILABLE ROUTINES FOR BETTER
C     PERFORMANCE, SUCH AS 'RMOVE' HAVE BEEN COLLECTED IN THIS ROUTINE.
C
C*****	    ENTRY RMOVE(X,Y,N)
C*****	RMOVE was changed to a separate subroutine.
C
C     MOVE 'N' COMPONENTS OF ARRAY 'Y' TO ARRAY 'X'.
C     THIS LOGIC WILL HANDLE OVERLAP CASES ONLY FOR X(M)=X(M+N1)
C     WITH N1>0, M=1,N.
C*****      REAL X(N),Y(N)
C
C      DO 10 I=1,N
C10    X(I)=Y(I)
C     RETURN
C
C     ENTRY IMOVE(IX,IY,N)
C*****	    ENTRY IMOVE(X,Y,N)
C***** IMOVE was changed to a separate subroutine.
C
C     SPECIAL ENTRY FOR INTEGER ARRAYS FOR MACHINES HAVING DIFFERENT
C     WORD SIZES FOR INTEGERS THAN REALS.
C     INTEGER IX(N),IY(N)
C
C      IF(N.LE.0) RETURN
C     USE CYBER BUILT-IN ROUTINE
C*****      CALL MOVLEV(Y,X,N)
C
C     THE EQUIVALENT FORTRAN CODING FOR THIS ROUTINE ARE THE FOLLOWING
C     TWO LINES.
C      DO 20 I=1,N
C20    IX(I)=IY(I)
C
C*****	RETURN
C*****	ENTRY SYSSEC(SECS)
C
C     OBTAIN CURENT TIME OF DAY IN SECONDS FROM SYSTEM ROUTINE
C      SECS=0
C     SYSTEM CALL NEEDED FOR THE CYBER NOS OP SYS
C*****      CALL SECOND(SECS)
C     SYSTEM CALL NEEDED FOR CRAY COMPUTER UNDER CTSS
C      CALL TIMEUSED(ICPU,IO,ISYS,MEM)
C      SECS=ICPU/1.E6
C*****	RETURN
C
C*****	ENTRY SYSTIM(TIM)
C
C     OBTAIN CURENT TIME OF DAY AS A CHARACTER STRING FROM SYSTEM ROUTI
C      TIM='TIME'
C     SYSTEM CALL NEEDED FOR THE CYBER NOS OP SYS
C*****      CALL TIME(TIM)
C     SYSTEM CALL NEEDED FOR THE CRAY CTSS OP SYS
C     CALL TIMEDATE(TIM,DATX,MACH)
C*****	RETURN
C
C*****	ENTRY SYSDAT(DAT)
C
C     OBTAIN CURENT DATE AS A CHARACTER STRING FROM SYSTEM ROUTINE
C      DAT='DATE'
C     SYSTEM CALL NEEDED FOR THE CYBER NOS OP SYS
C*****      CALL DATE(DAT)
C     SYSTEM CALL NEEDED FOR THE CRAY CTSS OP SYS
C     CALL TIMEDATE(TIMX,DAT,MACH)
C*****	RETURN
C     STANDARD SYSTEM INITIALIZATION TO BE CALLED FROM MAIN PROGRAM
      ENTRY SYSINI
C     SYSTEM INITIALIZATION CALL FOR CTSS ON CRAY COMPUTER
C     CALL LINK("READ5,UNIT5=TTY,PRINT6,UNIT6=TTY//")
      RETURN
C
C*****	ENTRY MORTEM
C*****	STOP 'POST MORTEM'
C ----------------------------------------------------------------------!-------
      END SUBROUTINE DFSYS_test 
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      FUNCTION ND2B_test (INDEC) ! 27 Feb 2006
C ----------------------------------------------------------------------!-------
C
C     CONVERT A NUMERIC DECIMAL DIGIT STRING TO A BIT STRING
C     WITH EACH BIT SET BY A CORRESPONDING DECIMAL DIGIT.
C
      ND2B=0
      IF(INDEC.LE.0) RETURN
      M=INDEC
      MB=1
C
10    IF(MOD(M,10).NE.0) ND2B=ND2B+MB
      MB=MB*2
      M=M/10
      IF(M.GT.0) GO TO 10
C ----------------------------------------------------------------------!-------
      END FUNCTION ND2B_test 
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      FUNCTION GET1_test(Z)
C ----------------------------------------------------------------------!-------
C     FUNCTIONALLY THE SAME AS 'GET' PROGRAM, SEE DOCUMENTATION THERE
C     NEEDED BECAUSE RECEIVER MODELS WILL CALL GET TO OBTAIN TERRAIN VA
C     BUT THEY THEMSELVES ARE CALLED VIA GET. HENCE HAVE RE-ENTRANCE
C     PROBLEM.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      double precision d1mach                                     ! 10 Nov 2009
      DOUBLE PRECISION PF(10),PG(10),z(12)
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 "RKAM" INSERTED HERE
      DOUBLE PRECISION KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)
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     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 "RMACH" INSERTED HERE
      COMMON/CRMACH/rmach(9)
      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/CGET/ZERO
C
      IENTRY=1
      GO TO 5
C
C ----------------------------------------------------------------------!-------
      ENTRY DOT1_test(Z)
C ----------------------------------------------------------------------!-------
      IENTRY=2
5     CONTINUE
C 26 Aug 2008      IF(ZERO.EQ.0.d0) ZERO=EARTHR*RMACH(3)*2.d0
      IF(ZERO.EQ.0.d0) ZERO=EARTHR*D1MACH(3)*2.d0                  ! 26 Aug 2008
C
      IF(Z(NTIME).EQ.RKTIME) GO TO 10
C
      IF(Z(NSELECT).EQ.USELECT) CALL SURFACE
      IF(Z(NSELECT).EQ.FSELECT) CALL RECEVER
      IF(Z(NSELECT).EQ.GSELECT) CALL TOPOG
      Z(NTIME)=RKTIME
C     REMOVE MACHINE ROUND OFF NOISE FROM EXACT RECEIVER LOCATIONS
      IF(dABS(Z(NZ)).LE.ZERO) Z(NZ)=0.d0
C
10    IF(IENTRY.NE.1) GO TO 20
      GET1=Z(NZ)
      RETURN
C
20    GET1=Z(NPZR)*DRDT(1)+Z(NPZTH)*DRDT(2)+Z(NPZPH)*DRDT(3)
      dot1=get1
      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION GET1_test
C ----------------------------------------------------------------------!-------
C
C
