C file S5.FOR
C ----------------------------------------------------------------------!-------
! 25 Aug 2011 - Changes recommended by Larry Baker
! 23 Dec 2010 - Added new subroutine FTANH3, that is like FTANH, but also 
!               furnishes derivatives of T with respect to profile height values
! 22 Jun 2010 - Added comments
! 26 Apr 2010 - Added SUBROUTINES QUARTIC and CUBIC
C 16 Mar 2010 - Increased dimension of ..GP.. arrays to 611
C 15 Mar 2010 - Added SUBROUTINE PLTNOTE
C 15 Mar 2010 - Added common blocks B21,B22,B23,B24 to READW
C 10 Nov 2009 - Added:     double precision d1mach
C 30 Apr 2009 - replaced REAL*8 with DOUBLE PRECISION
C 20 Mar 2009 - Corrected units conversion of a frequency as a period in seconds
C 19 Mar 2009 - Added units conversion for a central Earth angle in nautical
C               miles
C 10 Sep 2008 - Made data types consistent.
C 27 Aug 2008 - Changed mach to mach_test
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 26 Aug 2008 - Added I1mach.for, R1mach.for, and D1mach.for
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 12 Mar 2008 - Correction to common block /HDRC/
C 11 Mar 2008 - Added some comments to DFSYS
C  7 Mar 2008 - Updated the machine constants in rmach()
C  5 Mar 2008 - Recode Subroutines SYSSEC, SYSTIM, and SYSDAT to use
C                standard Fortran 90 instrinsic Subroutine DATE_AND_TIME
C                instead of Microsoft's GETTIM and GETDAT, courtesy of:
C                Larry Baker, US Geological Survey, 650-329-5608, baker@usgs.gov
C 23 Feb 2007 - Changed input units conversion for amplitude in dB
C  6 Sep 2006 - Added largest floating-point number to RMACH
C 14 Jul 2006 - Added debug print statements, now commented out
C 13 Jul 2006 - Added debug print statements, now commented out
C 24 Mar 2006 - Moved GET from contour1
C 24 Mar 2006 - Moved GAMANG from S4
C 24 Mar 2006 - Moved GET, WCHANGE, and RENORM from S2
C 6 Mar 2006 - Gave values to setw and setovr
C 1 Mar 2006 - Brought in some subprograms from S0, S1, S2, and S3
C ----------------------------------------------------------------------!-------
Contains DFSYS from S0
Contains I1mach, R1mach, D1mach
Contains READW1, GTUNIT, SREAD1, READW from S0
Contains CLEAR, ND2B, UCON from S1
Contains GET1, GET, ITEST, ITOC, SETTRC, WCHANGE, RENORM, SET2, ISET2 from S2
Contains PUTDES, NUMSTG, SFILL, STRIM, RERR, RERROR, STOPIT, PUTKST, PUTKCT, 
C OPNREP, ZAPFIL, OVERRD, SFILTR, ALCOSH, GAUSEL, FTANH, FTANH2 from S3, FTANH3
Contains GAMANG from S4
Contains PLTNOTE, QUARTIC, CUBIC
C ----------------------------------------------------------------------!-------
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 10 Jun 2003 (not marked) changed lines 381 (SUBROUTINE DFSYS)
C 10 Jun 2003 (not marked) changed lines 927, 928, 935 (FUNCTION READW)
C ----------------------------------------------------------------------!-------
C 2 Jun 2004 - fixed two bugs related to character string length
C ----------------------------------------------------------------------!-------
C 10 Jun 2003 - line 837 (SUBROUTINE SET2) ENTRY ISET2 to SUB ISET2 & rewrote
C ----------------------------------------------------------------------!-------
C 3 January 2006. Added SUBROUTINE FTANH2 to support 2nd derivatives.
C 17 Sep 2004.Changed tanh to dtanh in SUBROUTINE FTANH
C  5 Jun 2004.  Separated SUBROUTINE OCNHD into its own file.
C  2 Jun 2004.  Added printing of pressure and molecular weight to 
C               ENTRY HEADER1 in SUBROUTINE OCNHD
C 28 April 2004. Modified SUBROUTINE FTANH.
C 25 April 2004. Modified SUBROUTINE FTANH.
C 21 April 2004. Added SUBROUTINE FTANH at end of file.
C Sometime in 2003, changed line 17 of FUNCTION RERR (ERT = ...).
C 6 May 2003. Added 3 lines at line 58 in ENTRY HEADER1.
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
C from S0
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE DFSYS
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  11 Mar 2008.  The main idea here is that these routines may depend on
c  the local installation or system, so the user may have to substitute 
c  local versions of some of these routines.  The versions here are those
c  we use in our local system.  Some of these may work in other systems,
c  also.
 
      implicit double precision (A-H,O-Z)
      CHARACTER DAT*(*),DATX*10
      CHARACTER TIM*(*),TIMX*10
      character pc*1,typ*20
      common/crmach/rmach(9)
      integer ix(*),iy(*)
      double precision x(*),y(*)
      double precision d1mach                                     ! 10 Nov 2009

C    Declaration of new local variables:                           !  5 Mar 2008

      CHARACTER ZONX*10                                            !  5 Mar 2008
      INTEGER VALUES(8)                                            !  5 Mar 2008
      DOUBLE PRECISION MSEC                                        !  5 Mar 2008

C ----------------------------------------------------------------------!-------
      ENTRY RMOVE(X,Y,N)
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     replaces CYBER MOVLEV routine
      DO 10 I=1,N
10    X(I)=Y(I)
      RETURN
C ----------------------------------------------------------------------!-------
      ENTRY IMOVE(IX,IY,N)
C ----------------------------------------------------------------------!-------
C     SPECIAL ENTRY FOR INTEGER ARRAYS FOR MACHINES HAVING DIFFERENT
C     WORD SIZES FOR INTEGERS THAN REALS.
C     replaces CYBER BUILT-IN ROUTINE:
C     CALL MOVLEV(Y,X,N)
       DO 20 I=1,N
 20    IX(I)=IY(I)
      RETURN
 
 
C ----------------------------------------------------------------------!-------
      ENTRY SYSSEC(SECS)
C ----------------------------------------------------------------------!-------
c      ms fortran call to get system time for run times
C     call gettim(ihr,imin,isec,i100th)
C     xhr=float(ihr)
C     xmin=float(imin)
C     secs=float(isec)
C     x100th=float(i100th)/100.d0
C     secs=secs+x100th

C     Standard Fortran 90 intrinsic Subroutine DATE_AND_TIME        ! 5 Mar 2008
      Call DATE_AND_TIME( DATX, TIMX, ZONX, VALUES )                ! 5 Mar 2008
C     Date is returned as 'CCYYMMDD'                                ! 5 Mar 2008
C     Time is returned as 'hhmmss.sss'                              ! 5 Mar 2008
C     Zone is returned as '+hhmm' or '-hhmm'                        ! 5 Mar 2008
C     Values are:                                                   ! 5 Mar 2008
C        1: Year, 2: Month of year, 3: Day of month, 4: Minutes UTC difference,
C        5: Hour of the day, 6: Minutes of the hour, 7: Seconds of the minute,
C        8: Milliseconds of the second                              ! 5 Mar 2008
      xhr=float(VALUES(5))                                          ! 5 Mar 2008
      xmin=float(VALUES(6))                                         ! 5 Mar 2008
      secs=float(VALUES(7))                                         ! 5 Mar 2008
      msec=float(VALUES(8))/1000.d0                                 ! 5 Mar 2008
      secs=secs+msec                                                ! 5 Mar 2008
      secs=secs+60.d0*xmin
      secs=secs+3600.d0*xhr
      RETURN

 
C ----------------------------------------------------------------------!-------
      ENTRY SYSTIM(TIM)
C ----------------------------------------------------------------------!-------
c     ms fortran system time call for time labels
C     call gettim(ihr,imin,iday,i100th)
C     write (tim,'(2x,i2,1hh,i2,1hm,2x)') ihr,imin

C     Standard Fortran 90 intrinsic Subroutine DATE_AND_TIME        ! 5 Mar 2008
      Call DATE_AND_TIME( DATX, TIMX )                              ! 5 Mar 2008
C     Date is returned as 'CCYYMMDD'                                ! 5 Mar 2008
C     Time is returned as 'hhmmss.sss'                              ! 5 Mar 2008
      TIM = '  ' // TIMX(1:2) // 'h' // TIMX(3:4) // 'm'            ! 5 Mar 2008
      RETURN

C ----------------------------------------------------------------------!-------
      ENTRY SYSDAT(DAT)
C ----------------------------------------------------------------------!-------
c     ms fortran call to get system date for labels
C     call getdat(iyr,imon,iday)
C     write (dat,'(i2,1h/,i2,1h/,i4)') imon,iday,iyr

C     Standard Fortran 90 intrinsic Subroutine DATE_AND_TIME        ! 5 Mar 2008
      Call DATE_AND_TIME( DATX, TIMX )                              ! 5 Mar 2008
C     Date is returned as 'CCYYMMDD'                                ! 5 Mar 2008
C     Time is returned as 'hhmmss.sss'                              ! 5 Mar 2008
      DAT = DATX(5:6) // '/' // DATX(7:8) // '/' // DATX(1:4)       ! 5 Mar 2008
      RETURN

C ----------------------------------------------------------------------!-------
      ENTRY MORTEM
C ----------------------------------------------------------------------!-------
c     can add system post-mortem dump here
      STOP 'POST MORTEM'
c
c
C ----------------------------------------------------------------------!-------
      entry pcodes (pc,typ)
C ----------------------------------------------------------------------!-------
c     insert printer codes for form feed (pc), typography, etc.
c     the ones defined here are for the HP LaserJet 2
c     using line printer font with 8 lpi. 
c     LaserJet 1 uses A or B cartridge and prints in landscape mode.
c     pc=^L will work for most printers
      pc='16#0C'
! 25 Aug 2011      write (typ,'(a1,a4,a1,a8)')16#1B,'&l8D',16#1B,'(s16.66H'
      write (typ,'(a1,a4,a1,a8)')CHAR(27),'&l8D',CHAR(27),'(s16.66H'!25 Aug 2011
      return
c
c
C ----------------------------------------------------------------------!-------
      entry mach_test
C ----------------------------------------------------------------------!-------
C                                                                       DFCN0130
C     SINGLE PRECISION MACHINE CONSTANTS                                DFCN0140
C      (SEE DIGITAL SIGNAL PROCESSING, IEEE PRESS 1979 P S-7)           DFCN0150
C                                                                       DFCN0160
C     THE CONSTANTS ARE DERIVED FROM THE FOLLOWING                      DFCN0170
C     FORM FOR FLOATING POINT NUMBERS:                                  DFCN0180
C                                                                       DFCN0190
C        SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )                   DFCN0200
C     FOR MOST MACHINES THE BASE B=2.                                   DFCN0210
C     FOR THE CYBER 60-BIT MACHINES                                     DFCN0220
C        T=48, EMIN=-974 AND EMAX=1070.                                 DFCN0230
C     FOR THE PDP-11 AND VAX-11 MACHINES                                DFCN0240
C        T=24, EMIN=-127 AND EMAX=127                                   DFCN0250
C                                                                       DFCN0260
C     RMACH(1)  SMALLEST POSITIVE MAGNITUDE = B**(EMIN-1)               DFCN0270
C     RMACH(2)  LARGEST MAGNITUDE           = B**EMAX*(1-B**(-T))       DFCN0280
C     RMACH(3)  SMALLEST RELATIVE SPACING   = B**(-T)                   DFCN0290
C     RMACH(4)  LARGEST RELATIVE SPACING    = B**(1-T)                  DFCN0300
C     RMACH(5)                              = LOG10(B=2)                DFCN0310
C                                                                       DFCN0320
c     replaces the DFCNST subroutine,
c     since these are the only machine constants used by HARPO
C Double-precision constants: ! 7 Mar 2008
C     B=2 ! 7 Mar 2008
C     EMAX =  1023 ! 7 Mar 2008 ! EMAX =  1024                     ! 26 Aug 2008
C     EMIN = -1022 ! 7 Mar 2008 ! EMIN = -1021                     ! 26 Aug 2008
C 26 Aug 2008      rmach(1)=2.d0**(-1022)!smalleest positive magnitude,7 Mar 2008
C 26 Aug 2008      rmach(2)=(2.d0**(1023.d0)*(1.d0-2.d0**(-53.d0)))*2.d0 ! largest 
C                                              floating-point number, 6 Sep 2006
C 26 Aug 2008      rmach(3)=2.d0**(-53)  ! smallest relative spacing
C 26 Aug 2008      rmach(4)=2.d0**(-52)  ! largest relative spacing,  7 Mar 2008
C 26 Aug 2008      rmach(5)=dlog10(2.d0) 
C 26 Aug 2008      rmach(6)=2147483647.d0  ! largest 4-byte integer
      rmach(1)=d1mach(1)                                           ! 26 Aug 2008
      rmach(2)=d1mach(2)                                           ! 26 Aug 2008
      rmach(3)=d1mach(3)                                           ! 26 Aug 2008
      rmach(4)=d1mach(4)                                           ! 26 Aug 2008
      rmach(5)=d1mach(5)                                           ! 26 Aug 2008
      rmach(6)=I1mach(9)                                           ! 26 Aug 2008
      return
c
C ----------------------------------------------------------------------!-------
      END SUBROUTINE DFSYS
C ----------------------------------------------------------------------!-------
 
 
 
C ----------------------------------------------------------------------!-------
C later additions
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
C 26 Aug 2008 - Added I1mach.for, R1mach.for, D1mach.for
C ----------------------------------------------------------------------!-------
!DECK I1MACH
C ----------------------------------------------------------------------!-------
      INTEGER FUNCTION I1MACH (I)
      IMPLICIT NONE
      INTEGER :: I
      REAL :: X
      DOUBLE PRECISION :: XX
!***BEGIN PROLOGUE  I1MACH
!***PURPOSE  Return integer machine dependent constants.
!***LIBRARY   SLATEC
!***CATEGORY  R1
!***TYPE      INTEGER (I1MACH-I)
!***KEYWORDS  MACHINE CONSTANTS
!***AUTHOR  Fox, P. A., (Bell Labs)
!           Hall, A. D., (Bell Labs)
!           Schryer, N. L., (Bell Labs)
!***DESCRIPTION
!
!   I1MACH can be used to obtain machine-dependent parameters for the
!   local machine environment.  It is a function subprogram with one
!   (input) argument and can be referenced as follows:
!
!        K = I1MACH(I)
!
!   where I=1,...,16.  The (output) value of K above is determined by
!   the (input) value of I.  The results for various values of I are
!   discussed below.
!
!   I/O unit numbers:
!     I1MACH( 1) = the standard input unit.
!     I1MACH( 2) = the standard output unit.
!     I1MACH( 3) = the standard punch unit.
!     I1MACH( 4) = the standard error message unit.
!
!   Words:
!     I1MACH( 5) = the number of bits per integer storage unit.
!     I1MACH( 6) = the number of characters per integer storage unit.
!
!   Integers:
!     assume integers are represented in the S-digit, base-A form
!
!                sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
!
!                where 0 .LE. X(I) .LT. A for I=0,...,S-1.
!     I1MACH( 7) = A, the base.
!     I1MACH( 8) = S, the number of base-A digits.
!     I1MACH( 9) = A**S - 1, the largest magnitude.
!
!   Floating-Point Numbers:
!     Assume floating-point numbers are represented in the T-digit,
!     base-B form
!                sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
!
!                where 0 .LE. X(I) .LT. B for I=1,...,T,
!                0 .LT. X(1), and EMIN .LE. E .LE. EMAX.
!     I1MACH(10) = B, the base.
!
!   Single-Precision:
!     I1MACH(11) = T, the number of base-B digits.
!     I1MACH(12) = EMIN, the smallest exponent E.
!     I1MACH(13) = EMAX, the largest exponent E.
!
!   Double-Precision:
!     I1MACH(14) = T, the number of base-B digits.
!     I1MACH(15) = EMIN, the smallest exponent E.
!     I1MACH(16) = EMAX, the largest exponent E.
!
!   To alter this function for a particular environment, the desired
!   set of DATA statements should be activated by removing the C from
!   column 1.  Also, the values of I1MACH(1) - I1MACH(4) should be
!   checked for consistency with the local operating system.
!
!***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
!                 a portable library, ACM Transactions on Mathematical
!                 Software 4, 2 (June 1978), pp. 177-188.
!***ROUTINES CALLED  (NONE)
!***REVISION HISTORY  (YYMMDD)
!   750101  DATE WRITTEN
!   960411  Modified for Fortran 90 (BE after suggestions by EHG).   
!   980727  Modified value of I1MACH(6) (BE after suggestion by EHG).   
!***END PROLOGUE  I1MACH
!
      X  = 1.0      
      XX = 1.0D0

      SELECT CASE (I)
        CASE (1)
          I1MACH = 5 ! Input unit
        CASE (2)
          I1MACH = 6 ! Output unit
        CASE (3)
          I1MACH = 0 ! Punch unit is no longer used
        CASE (4)
          I1MACH = 0 ! Error message unit
        CASE (5)
          I1MACH = BIT_SIZE(I)
        CASE (6)
          I1MACH = 4            ! Characters per integer is hopefully no
                                ! longer used. 
                                ! If it is used it has to be set manually.
                                ! The value 4 is correct on IEEE-machines.
        CASE (7)
          I1MACH = RADIX(1)
        CASE (8)
          I1MACH = BIT_SIZE(I) - 1
        CASE (9)
          I1MACH = HUGE(1)
        CASE (10)
          I1MACH = RADIX(X)
        CASE (11)
          I1MACH = DIGITS(X)
        CASE (12)
          I1MACH = MINEXPONENT(X)
        CASE (13)
          I1MACH = MAXEXPONENT(X)
        CASE (14)
          I1MACH = DIGITS(XX)
        CASE (15)
          I1MACH = MINEXPONENT(XX)
        CASE (16)
          I1MACH = MAXEXPONENT(XX) 
        CASE DEFAULT
          WRITE (*, FMT = 9000)
 9000     FORMAT ('1ERROR    1 IN I1MACH - I OUT OF BOUNDS')
          STOP
        END SELECT
      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION I1MACH 
C ----------------------------------------------------------------------!-------
 
 
 
 
C ----------------------------------------------------------------------!-------
C 26 Aug 2008 - Added I1mach.for, R1mach.for, D1mach.for
C ----------------------------------------------------------------------!-------
!DECK R1MACH
C ----------------------------------------------------------------------!-------
      REAL FUNCTION R1MACH (I)
      IMPLICIT NONE
      INTEGER :: I
      REAL :: B, X
!***BEGIN PROLOGUE  R1MACH
!***PURPOSE  Return floating point machine dependent constants.
!***LIBRARY   SLATEC
!***CATEGORY  R1
!***TYPE      SINGLE PRECISION (R1MACH-S, D1MACH-D)
!***KEYWORDS  MACHINE CONSTANTS
!***AUTHOR  Fox, P. A., (Bell Labs)
!           Hall, A. D., (Bell Labs)
!           Schryer, N. L., (Bell Labs)
!***DESCRIPTION
!
!   R1MACH can be used to obtain machine-dependent parameters for the
!   local machine environment.  It is a function subprogram with one
!   (input) argument, and can be referenced as follows:
!
!        A = R1MACH(I)
!
!   where I=1,...,5.  The (output) value of A above is determined by
!   the (input) value of I.  The results for various values of I are
!   discussed below.
!
!   R1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
!   R1MACH(2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
!   R1MACH(3) = B**(-T), the smallest relative spacing.
!   R1MACH(4) = B**(1-T), the largest relative spacing.
!   R1MACH(5) = LOG10(B)
!
!   Assume single precision numbers are represented in the T-digit,
!   base-B form
!
!              sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
!
!   where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), and
!   EMIN .LE. E .LE. EMAX.
!
!   The values of B, T, EMIN and EMAX are provided in I1MACH as
!   follows:
!   I1MACH(10) = B, the base.
!   I1MACH(11) = T, the number of base-B digits.
!   I1MACH(12) = EMIN, the smallest exponent E.
!   I1MACH(13) = EMAX, the largest exponent E.
!
!
!***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
!                 a portable library, ACM Transactions on Mathematical
!                 Software 4, 2 (June 1978), pp. 177-188.
!***ROUTINES CALLED  XERMSG
!***REVISION HISTORY  (YYMMDD)
!   790101  DATE WRITTEN
!   960329  Modified for Fortran 90 (BE after suggestions by EG)      
!***END PROLOGUE  R1MACH
!      
      X = 1.0
      B = RADIX(X)
      SELECT CASE (I)
        CASE (1)
          R1MACH = B**(MINEXPONENT(X)-1) ! the smallest positive magnitude.
        CASE (2)
          R1MACH = HUGE(X)               ! the largest magnitude.
        CASE (3)
          R1MACH = B**(-DIGITS(X))       ! the smallest relative spacing.
        CASE (4)
          R1MACH = B**(1-DIGITS(X))      ! the largest relative spacing.
        CASE (5)
          R1MACH = LOG10(B)
        CASE DEFAULT
          WRITE (*, FMT = 9000)
 9000     FORMAT ('1ERROR    1 IN R1MACH - I OUT OF BOUNDS')
          STOP
      END SELECT
      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION R1MACH 
C ----------------------------------------------------------------------!-------
 
 
 
 
C ----------------------------------------------------------------------!-------
C 26 Aug 2008 - Added I1mach.for, R1mach.for, D1mach.for
C ----------------------------------------------------------------------!-------
!DECK D1MACH
C ----------------------------------------------------------------------!-------
      DOUBLE PRECISION FUNCTION D1MACH (I)
      IMPLICIT NONE
      INTEGER :: I
      DOUBLE PRECISION :: B, X
!***BEGIN PROLOGUE  D1MACH
!***PURPOSE  Return floating point machine dependent constants.
!***LIBRARY   SLATEC
!***CATEGORY  R1
!***TYPE      SINGLE PRECISION (D1MACH-S, D1MACH-D)
!***KEYWORDS  MACHINE CONSTANTS
!***AUTHOR  Fox, P. A., (Bell Labs)
!           Hall, A. D., (Bell Labs)
!           Schryer, N. L., (Bell Labs)
!***DESCRIPTION
!
!   D1MACH can be used to obtain machine-dependent parameters for the
!   local machine environment.  It is a function subprogram with one
!   (input) argument, and can be referenced as follows:
!
!        A = D1MACH(I)
!
!   where I=1,...,5.  The (output) value of A above is determined by
!   the (input) value of I.  The results for various values of I are
!   discussed below.
!
!   D1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
!   D1MACH(2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
!   D1MACH(3) = B**(-T), the smallest relative spacing.
!   D1MACH(4) = B**(1-T), the largest relative spacing.
!   D1MACH(5) = LOG10(B)
!
!   Assume single precision numbers are represented in the T-digit,
!   base-B form
!
!              sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
!
!   where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), and
!   EMIN .LE. E .LE. EMAX.
!
!   The values of B, T, EMIN and EMAX are provided in I1MACH as
!   follows:
!   I1MACH(10) = B, the base.
!   I1MACH(11) = T, the number of base-B digits.
!   I1MACH(12) = EMIN, the smallest exponent E.
!   I1MACH(13) = EMAX, the largest exponent E.
!
!
!***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
!                 a portable library, ACM Transactions on Mathematical
!                 Software 4, 2 (June 1978), pp. 177-188.
!***ROUTINES CALLED  XERMSG
!***REVISION HISTORY  (YYMMDD)
!   790101  DATE WRITTEN
!   960329  Modified for Fortran 90 (BE after suggestions by EHG)      
!***END PROLOGUE  D1MACH
!      
      X = 1.0D0
      B = RADIX(X)
      SELECT CASE (I)
        CASE (1)
          D1MACH = B**(MINEXPONENT(X)-1) ! the smallest positive magnitude.
        CASE (2)
          D1MACH = HUGE(X)               ! the largest magnitude.
        CASE (3)
          D1MACH = B**(-DIGITS(X))       ! the smallest relative spacing.
        CASE (4)
          D1MACH = B**(1-DIGITS(X))      ! the largest relative spacing.
        CASE (5)
          D1MACH = LOG10(B)
        CASE DEFAULT
          WRITE (*, FMT = 9000)
 9000     FORMAT ('1ERROR    1 IN D1MACH - I OUT OF BOUNDS')
          STOP
      END SELECT
      RETURN
C ----------------------------------------------------------------------!-------
      END  FUNCTION D1MACH 
C ----------------------------------------------------------------------!-------


 
C ----------------------------------------------------------------------!-------
C from S0
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE READW1(AB,NWOK,MSET,MX,NTBL,ITBL,FRMTBL,GP)
C ----------------------------------------------------------------------!-------
C     HANDLES INPUT OF TABULAR DATA REQUIRED BY SOME MODELS.
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
C     COMMON DECK "CUCON" INSERTED HERE
      COMMON/UCONV/CNVV(5,4)
      CHARACTER PCV*3,CNVC*2
      COMMON/UCONC/PCV(5),CNVC(5,4)
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 "WWR" INSERTED HERE
      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))
      LOGICAL NWOK,AB,ucon,uucon
C
      CHARACTER LINE*80,ALPHA*3,FBUF*3
      PARAMETER (NVAR=6,ALPHA='  A',NACDE=-1,NFIELD=13)
      INTEGER IPV(NVAR),NTBL(10),ITBL(10),FRMTBL(10)
      CHARACTER*8 PU(NVAR)
      DOUBLE PRECISION GP(*),V(NVAR),CONV(NVAR)
C
C     THE DATA STRUCTURE EXPECTED BY THIS PROGRAM FOR A DATA BLOCK IS
C     SPECIFIED FOR A GIVEN MODEL BY THE THREE ARRAYS NTBL,ITBL,FRMTBL.
C     THE TYPE OF DATA EXPECTED FOR A GIVEN FORMAT GROUP IS TAKEN FROM
C     ARRAY FRMTBL. IT MUST HAVE AN ENTRY FOR A GIVEN TYPE OF FORMAT
C     WHICH IS EITHER ZERO, IF THAT FORMAT IS NOT ALLOWED OR AN ENTRY
C     WHOSE VALUE IS EQUAL TO THAT FORMAT TYPE. CURRENTLY ONLY TWO
C     FORMAT TYPES ARE ALLOWED, 1 FOR ALPHA AND 2 FOR NUMERIC.
C     TO ALLOW MORE FLEXIBLITY IN THE INPUT FILE THE NUMERIC FORMAT
C     TYPE IN TURN ALLOWS FOR 3 DIFFERENT DATA FORMATS, SEE BELOW.
C
C     OFFSETS INTO THE GENERAL ARRAY 'GP' FOR EACH FORMAT ARE GIVEN
C     IN THE ARRAY 'NTBL'. THE SPACING BETWEEN DATA VALUES IN EACH LINE
C     READ IS GIVEN IN ARRAY 'ITBL'. THIS SCHEME DOES NOT ALLOW FOR X,Y
C     ,Z CYCLES BUT RATHER SEPARATE ARRAYS OF X VALUES, THEN Y VALUES,
C     AND THEN Z VALUES(ETC. UP TO 'NVAR' VARIABLES).
C
      NWOK=.TRUE.
C
C     INITIALIZE ANY ALPHANUMERIC ARRAYS TO BLANKS
      dummy=1h 
      IF(FRMTBL(1).EQ.1) CALL SET2(GP(NTBL(1)),dummy,NTBL(2)-NTBL(1))
C
C     BEGIN MULTI-FORMAT LOOP
5     READ(NRYIND,'(A,A)',END=200) FBUF,LINE
      IF(NDEVTMP.GT.0) WRITE(NDEVTMP,'(A,A)') FBUF,LINE
c      print *, ' (zeroth) fbuf = ', fbuf ! 25 Mar 2009
c      print *, ' (zeroth) alpha = ', alpha ! 25 Mar 2009
      IF(FBUF.EQ.ALPHA) THEN
         NFRM=NACDE
      ELSE
         READ(FBUF,10) NFRM
c      print *, ' (zeroth) nfrm = ', nfrm ! 25 Mar 2009
      ENDIF
10    FORMAT(BZ,I3,A)
      IF(NFRM .EQ. 0) RETURN
C
C     FORMAT #1 IS ASSIGNED TO ALPHANUMERIC INPUT AND LARGER NUMBERS
C     FOR NUMERIC INPUT.
C
C     CHECK NOW FOR NUMERICAL FORMATS.
C     THESE ARE PROVIDED FOR EASE OF INPUT OF TABULATED DATA
C     AND ARE PARTLY 'TRANSPARENT' TO THE MODEL INVOLVED. I.E.
C     THE USER HAS THE OPTION OF SPECIFING THE NUMBER OF DATA COLUMNS.
C     INPUT FORMAT NUMBER SPECIFIES THE NUMBER OF COLUMNS OF INPUT DATA
C     I.E. FOR VALUE 1 A SINGLE DATA COLUMN IS EXPECTED, FOR A VALUE
C     2 TWO DATA COLUMNS ARE EXPECTED, ETC.
C
C     DETERMINE FORMAT NUMBER
c      print *, ' (first) mset = ', mset ! 13 Jul 2006
c      print *, ' (first) nacde = ', nacde ! 13 Jul 2006
c      print *, ' (first) nfrm = ', nfrm ! 13 Jul 2006
c      print *, ' (first) mx = ', mx ! 13 Jul 2006
c      print *, ' (first) it = ', it ! 13 Jul 2006
      IF(NFRM.EQ.NACDE) THEN
         IT=1
C     ALPHANUMERIC FORMATS MUST BE FIRST IN LIST
C     SINCE NO OFFSET SPECIFIER IS ALLOWED IN A FULL ALPHA LINE
c      print *, ' FRMTBL(1) = ', FRMTBL(1)  ! 13 Jul 2006
         IF(FRMTBL(1).EQ.1) GO TO 100
      ELSE
         READ(LINE,'(I2)') IT
c      print *, ' (second) it = ', it ! 13 Jul 2006
         IT=MAX0(IT,2)
         IF(IT.LE.MX) GO TO 100
      ENDIF
C
c      print *, ' (second) nacde = ', nacde ! 13 Jul 2006
15    AB=.TRUE.
c      print *, ' (third) nacde = ', nacde ! 13 Jul 2006
c      print *, ' it = ', it ! 13 Jul 2006
      PRINT 20, MSET,NFRM,MX
20    FORMAT(' FOR SET',I5,' FORMAT NUMBER',I5,' EXCEEDS LIMIT OF',I5)
C      ERROR IN INPUT DATA , STOP HERE
      STOP
C
C     THE 'NTBL' TABLE PROVIDES A LIST OF OFFSETS TO VARIABLE GROUPS IN
C     MODEL ARRAYS, 'ITBL' GIVES ELEMENT SEPARATION TO ALLOW FOR 2-
C     DIMENSIONAL ARRAYS AND 'FRMTBL' SPECIFIES THE INPUT FORMAT.
C
100   N1=NTBL(IT)
      N3=ITBL(IT)
! 25 Aug 2011      N2=NTBL(IT+1)-N3
      N2=(N1+1)+N3*(((NTBL(IT+1)-(N1+1))/N3)-1)                    ! 25 Aug 2011
C
      IF(NFRM.NE.NACDE) GO TO 110
C
      READ(LINE,1010)((GP(I+J-1),I=N1,N2,N3),J=1,N3)
1010  FORMAT(10A8)
      GO TO 5
C
110   IF(NFRM.LT.1 .OR. NFRM.GT.NVAR) GO TO 200
C     HANDLE NUMERIC FORMATS
C
C     N1 IS STARTING ELEMENT WHICH WILL BE AFTER THE DATA COUNT
C     VALUE(NEED NELS*N3+1 TOTAL ELEMENTS)
      N1=N1+1
      N=0
C
C     # ELEMENTS IS EQUAL TO FORMAT NUMBER
      NELS=NFRM
C     READ UNITS CONVERSION SPECIFICATION LINE
      READ(LINE,1018) (PU(I),I=1,3)
      IF(PU(2).NE.' ') READ(LINE,2024) DELIM
      IF(PU(2).EQ.' ') READ(LINE,1020) DELIM,(PU(I),I=1,NELS)
C
C  ALLOW FOR BLANK CONVERSION LINE PRODUCED BY EIGENRAY ROUTINE
      IF(DELIM.EQ.0) READ(LINE,2024) DELIM
C
C     READ NEXT LINE AND EXTRACT UNIT CONVERSION SPECS IF PRESENT
C     IN FLOATING FORMAT
      READ(NRYIND,'(A)') LINE
      CALL GTUNIT(LINE,PU,NFIELD,NELS, NCARY)
C     CARRY COUNT IS ZERO IF LINE CONTAINED UNITS SPECIFIERS
      IF(NCARY.EQ.0.AND.NDEVTMP.GT.0) WRITE(NDEVTMP,'(A)') LINE
C
1018  FORMAT(10(A,5X))
1020  FORMAT(BZ,3X,G10.3,8X,10A)
2024  FORMAT(T43,G10.3)
1025  FORMAT(BZ,10G13.6)
C
C     LOOK UP CONVERSION CONSTANTS
      DO 2030 I=1,NELS
2030  uucon= UCON(IPV(I),PU(I),CONV(I))
C
C     DELIMITOR TESTING REQUIRES THE INNER LOOP BE EXECUTED
C     AT LEAST ONCE
      NX=MAX0(1,N3)
      N2=MAX0(N2,N1)
C     BEGIN SINGLE DATA FORMAT LOOP
      DO 150 J=0,10000
         IF(NELS.GT.1) THEN
            IF(NCARY.EQ.0) READ(NRYIND,'(A)') LINE
            NCARY=0
            IF(NDEVTMP.GT.0) WRITE(NDEVTMP,'(A)') LINE
            READ(LINE,1025) (V(I),I=1,NELS)
         ENDIF
         K=0
         DO 145 I=N1+J,N2+J,NX
               K=K+1
               IF(NELS.EQ.1) THEN
                  IF(NCARY.EQ.0) READ(NRYIND,'(A)',END=10000) LINE
                  NCARY=0
                  IF(NDEVTMP.GT.0) WRITE(NDEVTMP,'(A)') LINE
C     FOR SINGLE COLUMN CASE READ NOW
                  READ(LINE,1020) V(K)
               ENDIF
               IF(V(K) .EQ. DELIM) GO TO 160
C
C           APPLY ANY SCALING CONVERSIONS
               IF(CONV(K).GE.0.d0) THEN
                  V(K)=V(K)*CONV(K)
               ELSEIF(V(K).NE.0.d0) THEN
                  V(K)=-CONV(K)/V(K)
      ENDIF
               IF(I.GE.NTBL(IT+1)) THEN                            ! 25 Aug 2011
                  WRITE (*,*) 'GP TABLE OVERFLOW IN DATA SET ', MSET,
     1                        ': NPT=', N+1, ', ROW=', J+1, ', COL=', K
                  CALL RERROR('READW1','OVERFLOW',FLOAT(I))
               ENDIF                                               ! 25 Aug 2011
               IF(J.LT.N3) GP(I)=V(K)
145            N=N+1
150   CONTINUE
C
10000 IF(V(1).NE.DELIM) CALL RERROR('READW1','NO DELIM',FLOAT(MSET))
C
C     SAVE NUMBER OF VALUES AND ELEMENTS(AS FRACTIONAL PART)
160   IF(N1.GT.1) GP(N1-1)=N+NELS/10.

      WRITE(*,*) 'DATA SET ', MSET, ': NPTS=', N,                  ! 25 Aug 2011
     1           ', ROWS=', (N+NELS-1)/NELS, ', COLS=', NELS       ! 25 Aug 2011
C
C     CONTINUE MULTI-FORMAT LOOP
      GO TO 5
C
200   CALL RERROR('READW1','BAD FORMAT',FLOAT(MSET))
C
C ----------------------------------------------------------------------!-------
      END SUBROUTINE READW1
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE GTUNIT(LINE,PU,NFIELD,NVAR, LN)
C ----------------------------------------------------------------------!-------
C     GTUNIT IS USED TO EXTRACT UNIT CONVERSION SPECIFIERS FROM TABULAR
C     DATA INPUT.
C     CONVERSION SPECIFICATIONS ARE TAKEN FROM 'LINE' INPUT STRING.
C     A FLOATING FORMAT FALLING WITHIN SUCCESSIVE FIELDS OF
C     LENGTH 'NFIELD' IS PROVIDED FOR.
C     OUTPUT IS TO UNIT SPECIFIER ARRAY 'PU' FOR A MAXIMUM OF 'NVAR'
C     FIELDS. FOR A SUCCESS/FAIL RETURN VALUE, 'LN' GIVES THE NUMBER
C     OF NON-BLANK CHARACTERS USEFUL FOR UNIT DECODING, ZERO IMPLIES
C     NO USEFUL UNIT INFORMATION WAS FOUND.
C
 
      implicit double precision (a-h,o-z)
 
      PARAMETER (MAXLN=130)
      CHARACTER*(*) LINE,PU(NVAR),TMP*(MAXLN),TMP1*(MAXLN)
      INTEGER STRIM
C
      LN=STRIM(LINE)
      IF(LN.LT.1) THEN
         DO 5 I=1,NVAR
5           PU(I)=' '
         RETURN
      ENDIF
C
      TMP=' '
      CALL SFILTR(LINE,TMP,' 0123456789+-.E')
      IF(STRIM(TMP).LT.1) RETURN
C     SIGNIFY WE CONSUMED THIS LINE
      LN=0
C
      TMP1=LINE
C     PLACE INTER-SPECIFIER DELIMITORS
      DO 10 I=NFIELD,NFIELD*NVAR,NFIELD
         IF(TMP1(I:I).NE.' ') THEN
            IF(TMP1(I+1:I+1).NE.' ') THEN
C     IF NO ROOM FOR A BLANK, ITS AN ERROR
               CALL RERROR('GTUNIT','FORMAT ER1',FLOAT(I))
            ELSE
               TMP1(I+1:I+1)='#'
            ENDIF
         ELSE
            TMP1(I:I)='#'
         ENDIF
10    CONTINUE
C
C     NOW SQUEEZE OUT THE BLANKS
      CALL SFILTR(TMP1,TMP,' ')
C     HOPEFULLY SOMETHING LEFT
C
      NL=STRIM(TMP)
      ILIN=1
      DO 20 I=1,NVAR
C     CHECK FOR EMPTY CONVERSIONS
         IF(TMP(ILIN:ILIN).EQ.'#') THEN
            ILIN=ILIN+1
            PU(I)=' '
            GO TO 20
         ENDIF
C     GET NEXT UNIT TYPE
         PU(I)(1:2)=TMP(ILIN:ILIN+1)
C     GET NEXT PHYSICAL UNIT(BLANKS HAVE BEEN REMOVED)
C     SHORT PHRASE IF A BLANK WAS REMOVED
         IF(TMP(ILIN+3:ILIN+3).NE.'#') THEN
            IF(TMP(ILIN+4:ILIN+4).NE.'#') THEN
C     IF INCORRECT GROUPING OF LETTERS, ITS AN ERROR
               CALL RERROR('GTUNIT','FORMAT ER2',FLOAT(ILIN))
            ELSE
C     ADD TWO LETTER UNIT
               PU(I)(3:5)=' '//TMP(ILIN+2:ILIN+3)
               ILIN=ILIN+5
            ENDIF
         ELSE
C     ADD ONE LETTER UNIT
            PU(I)(3:5)='  '//TMP(ILIN+2:ILIN+2)
            ILIN=ILIN+4
         ENDIF
20    CONTINUE
C
C     SHOULD HAVE USED UP ALL THE TEXT BY NOW
      IF(ILIN.LE.NL) CALL RERROR('GTUNIT','FORMAT ER3',FLOAT(ILIN))
C ----------------------------------------------------------------------!-------
      END SUBROUTINE GTUNIT
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE SREAD1(AB,NWOK,NW)
C ----------------------------------------------------------------------!-------
 
      implicit double precision (a-h,o-z)
      logical ab,nwok
C     THIS IS THE FALLTHROUGH FEATURE FOR READW WHEN ENCOUNTERING
C     A NON STANDARD LABELED COMMON BLOCK. SREAD1 WILL ADD TABULAR
C     DATA TO THE DEFAULT GP ARRAY PGP AND INCLUDE INDEX OFFSET
C     VALUES TO A LOOKUP ARRAY 'LIST'.
C     COMMON DECK "CPROCFL" INSERTED HERE
      INTEGER PMX,PNTBL,PITBL,PFRMTBL,IDSP(10)
C     PARAMETER DECK "PGROUPS"
      PARAMETER (NCHPG1=11,NWPV=500,NSPGP=NCHPG1+2*NWPV+1)         ! 25 Aug 2011
      PARAMETER (MNGRP=9,MXGRP=69,MXLIST=MXGRP-MNGRP+2)
      COMMON/PROCFL/LIST(MXLIST)
      COMMON/PROCFL/PMX,PNTBL(10),PITBL(10),PFRMTBL(10),PGP(NSPGP)
      EQUIVALENCE (PGP,IDSP)
C
      CHARACTER ITOC*7
C
      INTEGER PPX    ,PQTBL(10),PLTBL(10),PIRMTBL(10)
C
      DATA PPX/2/
      DATA PQTBL/1,NCHPG1,NSPGP,7*0/
      DATA PLTBL/1,NWPV,8*0/
      DATA PIRMTBL/1,2,8*0/
C
      NWA=IABS(NW)-(MNGRP-2)
      IF(NWA.LT.2 .OR. NWA.GT.MXLIST) CALL STOPIT
     1  ( 'SREAD1 NW='//ITOC(NW)//',LIST(1)='//ITOC(LIST(1)) )
C
      CALL READW1(AB,NWOK,NW,PMX,PNTBL,PITBL,PFRMTBL,PGP)
C
C     UPDATE THE LIST, LIST(1)=MAXIMUM USED ELEMENT
      LIST(1)=MAX0(LIST(1),NWA)
      LIST(NWA)=PNTBL(1)
C
C     GET NEXT NUMERICAL INDEX FROM THE POINTER ARRAY
      NLN=PNTBL(2)
C
C     NUMBER OF VALUES IS INTEGRAL PART OF FIRST ELEMENT
      NPTS=PGP(NLN)
C     NUMBER OF ROWS IS KEPT AS FRACTIONAL PART
      NELS=(PGP(NLN)-NPTS)*10.+.5
! 25 Aug 2011      N2=NPTS/NELS
      N2=(NPTS+NELS-1)/NELS                                        ! 25 Aug 2011
      NMX=PITBL(2)
C     MAKE ADJUSTMENTS TO LENGTHS AND OFFSETS
! 25 Aug 2011      PITBL(2)=NMX-N2-(NCHPG1-1)
      NUP=NPTS+NCHPG1
      PNTBL(1)=PNTBL(1)+NUP
      PNTBL(2)=PNTBL(2)+NUP
      PITBL(2)=(PNTBL(3)-(PNTBL(2)+1))/PFRMTBL(2)                  ! 25 Aug 2011
      NLN=NLN+1
C
C     LOOP TO MOVE NELS-1 ROWS OF PGP(NMX,NELS) TO PGP(N2,NELS)
C     THUS MAKING PGP A CONTIGUOUS ARRAY
      NTO=NLN+N2
      NFRM=NLN+NMX
      DO 100 I=2,NELS
      IF(NTO.GT.0 .AND. NFRM+N2-1.LE.NSPGP)
     1        CALL RMOVE(PGP(NTO),PGP(NFRM),N2)
      NTO=NTO+N2
100   NFRM=NFRM+NMX
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY SETSRD
C ----------------------------------------------------------------------!-------
c%      CALL ISET2(LIST,0,MXLIST)
      do 66 i=1,mxlist
66    list(i)=0
      PMX=PPX
      CALL IMOVE(PNTBL,PQTBL,10)
      CALL IMOVE(PITBL,PLTBL,10)
      CALL IMOVE(PFRMTBL,PIRMTBL,10)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE SREAD1
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      FUNCTION   READW (XX)
C ----------------------------------------------------------------------!-------
C     READS W ARRAY MAKING ANY NECESSARY CONVERSION OF UNITS.
C     VARIABLE LENGTH 'TABULAR' DATA IS READ BY SUBROUTINE 'READW1'
C     WHICH IS CALLED FOR NEGATIVE W-ARRAY INDICES.
C
C     COMMON DECK "CONST" INSERTED HERE
 
      implicit double precision (a-h,o-z)
 
      COMMON/PCONST/CREF,RGAS,GAMMA
     1 ,MO2,MN2,MO,MN,MH2O,MO3                        ! 9 Mar 2007 ! 13 Mar 2008
      DOUBLE PRECISION MO2,MN2,MO,MN,MH2O,MO3                      ! 13 Mar 2008
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C     COMMON DECK "WWR" INSERTED HERE
      COMMON/WW/DUM(10),MAXW,W(1000)
      COMMON/IDENT/ID
      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 "B1" INSERTED HERE
      INTEGER UMX,UNTBL,UITBL,UFRMTBL,IDSU(20) 
      COMMON/B1/UMX,UNTBL(10),UITBL(10),UFRMTBL(10),UGP(611) 
      EQUIVALENCE (UGP,IDSU) 

C     COMMON DECK "B2" INSERTED HERE
      INTEGER DUMX,DUNTBL,DUITBL,DUFRMTB,IDSDU(20) 
      COMMON/B2/DUMX,DUNTBL(10),DUITBL(10),DUFRMTB(10),DUGP(611) 
      EQUIVALENCE (DUGP,IDSDU) 

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

C     COMMON DECK "B4" INSERTED HERE
      INTEGER DCMX,DCNTBL,DCITBL,DCFRMTB,IDSDC(20) 
      COMMON/B4/DCMX,DCNTBL(10),DCITBL(10),DCFRMTB(10),DCGP(611) 
      EQUIVALENCE (DCGP,IDSDC) 

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

C     COMMON DECK "B6" INSERTED HERE
      INTEGER DTMX,DTNTBL,DTITBL,DTFRMTB,IDSDT(20) 
      COMMON/B6/DTMX,DTNTBL(10),DTITBL(10),DTFRMTB(10),DTGP(611) 
      EQUIVALENCE (DTGP,IDSDT) 

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(611) 
      EQUIVALENCE (MGP,IDSM) 

C     COMMON DECK "B9" INSERTED HERE
      INTEGER GMX,GNTBL,GITBL,GFRMTBL,IDSG(20) 
      COMMON/B9/GMX,GNTBL(10),GITBL(10),GFRMTBL(10),GGP(611) 
      EQUIVALENCE (GGP,IDSG),(ANG,GGP(11)) 
C
C     COMMON DECK "B10" INSERTED HERE
      INTEGER DGMX,DGNTBL,DGITBL,DGFRMTB,IDSDG(20) 
      COMMON/B10/DGMX,DGNTBL(10),DGITBL(10),DGFRMTB(10),DGGP(611) 
      EQUIVALENCE (DGGP,IDSDG) 
C
C     COMMON DECK "B8" INSERTED HERE
      INTEGER RMX,RNTBL,RITBL,RFRMTBL,IDSR(20) 
      COMMON/B8/RMX,RNTBL(10),RITBL(10),RFRMTBL(10),RGP(611) 
      EQUIVALENCE (RGP,IDSR) 
C
C     COMMON DECK "CB11" INSERTED HERE
      INTEGER SMX,SNTBL,SITBL,SFRMTBL,IDSS(20) 
      COMMON/B11/SMX,SNTBL(10),SITBL(10),SFRMTBL(10),SGP(611) 
      EQUIVALENCE (SGP,IDSS),(ANS,SGP(11)) 

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

C     COMMON DECK "CB13" INSERTED HERE ! 31 Jan 2006
      INTEGER AMX,ANTBL,AITBL,AFRMTBL,IDSA(20) ! 31 Jan 2006
      COMMON/B13/AMX,ANTBL(10),AITBL(10),AFRMTBL(10),AGP(611) ! 31 Jan 2006
      EQUIVALENCE (AGP,IDSA),(ANA,AGP(11)) ! 31 Jan 2006

C     COMMON DECK "CB14" INSERTED HERE ! 31 Jan 2006
      INTEGER DAMX,DANTBL,DAITBL,DAFRMTB,IDSDA(20) ! 31 Jan 2006
      COMMON/B14/DAMX,DANTBL(10),DAITBL(10),DAFRMTB(10),DAGP(611) ! 31 Jan 2006
      EQUIVALENCE (DAGP,IDSDA),(ANDA,DAGP(11)) ! 31 Jan 2006

C     COMMON DECK "CB15" INSERTED HERE ! 31 Jan 2006
      INTEGER DMX,DNTBL,DITBL,DFRMTBL,IDSD(20) ! 31 Jan 2006
      COMMON/B15/DMX,DNTBL(10),DITBL(10),DFRMTBL(10),DGP(611) ! 31 Jan 2006
      EQUIVALENCE (DGP,IDSD),(AND,DGP(11)) ! 31 Jan 2006

C     COMMON DECK "CB16" INSERTED HERE ! 31 Jan 2006
      INTEGER DDMX,DDNTBL,DDITBL,DDFRMTB,IDSDD(20) ! 31 Jan 2006
      COMMON/B16/DDMX,DDNTBL(10),DDITBL(10),DDFRMTB(10),DDGP(611) ! 31 Jan 2006
      EQUIVALENCE (DDGP,IDSDD),(ANDD,DDGP(11)) ! 31 Jan 2006

C     COMMON DECK "CB17" INSERTED HERE
      INTEGER VMX,VNTBL,VITBL,VFRMTBL,IDSV(20) 
      COMMON/B17/VMX,VNTBL(10),VITBL(10),VFRMTBL(10),VGP(611) 
      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(611) 
      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(611) 
      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(611) 
      EQUIVALENCE (DPGP,IDSDP),(ANDP,DPGP(11)) 

C     COMMON DECK "CB21" INSERTED HERE                            ! 15 Mar 2010
      INTEGER MX21,NTBL21,ITBL21,FRMTB21,IDS21(20) 
      COMMON/B21/MX21,NTBL21(10),ITBL21(10),FRMTB21(10),GP21(611) 
      EQUIVALENCE (GP21,IDS21),(AN21,GP21(11)) 

C     COMMON DECK "CB22" INSERTED HERE                            ! 15 Mar 2010
      INTEGER MX22,NTBL22,ITBL22,FRMTB22,IDS22(20) 
      COMMON/B22/MX22,NTBL22(10),ITBL22(10),FRMTB22(10),GP22(611) 
      EQUIVALENCE (GP22,IDS22),(AN22,GP22(11)) 

C     COMMON DECK "CB23" INSERTED HERE                            ! 15 Mar 2010
      INTEGER MX23,NTBL23,ITBL23,FRMTB23,IDS23(20) 
      COMMON/B23/MX23,NTBL23(10),ITBL23(10),FRMTB23(10),GP23(611) 
      EQUIVALENCE (GP23,IDS23),(AN23,GP23(11)) 

C     COMMON DECK "CB24" INSERTED HERE                            ! 15 Mar 2010
      INTEGER MX24,NTBL24,ITBL24,FRMTB24,IDS24(20) 
      COMMON/B24/MX24,NTBL24(10),ITBL24(10),FRMTB24(10),GP24(611) 
      EQUIVALENCE (GP24,IDS24),(AN24,GP24(11)) 

      common/crmach/rmach(9)

C     COMMON DECK "HDR" INSERTED HERE
      CHARACTER*10 INITID*80,DAT                            ! ,TOD ! 12 Mar 2008
      common/hdr/sec,stsec                                 !Rev 7/13/90
      COMMON/HDRC/INITID,DAT,TOD                                   ! 12 Mar 2008

C 26 Aug 2008     1 ,(bigval,rmach(6))

      PARAMETER (MXCMTS=83)
      CHARACTER*38 WFRMT,WNOTES(2),bell*1,tod*10,face*1    !Rev 7/13/90
      character*2 numstg
      CHARACTER*80 LINEX
      CHARACTER*100 STMP1,STMP2
      INTEGER STRIM,WCOMTS(MXCMTS)
      CHARACTER*1 DEG,KM,NM,FEET,CYCLE,PER,MSKMS

C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVGRP,NDEVBIN
      LOGICAL NWOK,AB,uucon,ucon,setucon
C
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 ' CGRID2 ' inserted here		A. Newhall
	parameter (mx=7,my=7,mz=7)
	common/cgrid2/lxx7,nyy7,mzz7,xc7(7),yc7(7),zc7(7),cavg7(7)
	common/knots/zknot(11),yknot(11),xknot(11)
	common/coefs/cco(7),spleof(7,7),cbispl(7,7*7)
C
      CHARACTER PC*1,ID*80,typ*20
C
      INTEGER XMX,XNTBL(10),XITBL(10),XFRMTBL(10)
      DOUBLE PRECISION XGP(611)
C
      DATA XMX/2/
      DATA XNTBL/1,11,12,7*0/
      DATA XITBL/1,9*0/
      DATA XFRMTBL/1,2,8*0/
C
      DATA WCOMTS,WNOTES/MXCMTS*1, ' ',' INPUT OVERRIDDEN'/
*C
      AB=.FALSE.
      READW=0.d0
C
cc      print* , ' entering readw '
cc      print* , 'nryind = ', nryind
cc      print* , 'ndevtmp = ', ndevtmp
cc      print* , 'nfrmat = ', nfrmat
cc      print* , 'ndevgrp = ', ndevgrp
cc      print* , 'ndevbin = ', ndevbin
      IF(NDEVTMP.GT.0) REWIND NDEVTMP
C
      READ(NRYIND,'(a)',END=3) ID
      IF(NDEVTMP.GT.0) WRITE(NDEVTMP,'(a)') ID
c      print *,'id = ',id
C*      READ(NRYIND,1000,END=3) ID
C* 1000 FORMAT (BZ,10A8)
C*      IF(NDEVTMP.GT.0) WRITE(NDEVTMP,1000) ID
      GO TO 4
C
3     READW    =1.d0
c      print* , ' leaving readw because  NFRMAT = ', NFRMAT
      IF(NFRMAT.LT.0) RETURN
C
C%      PRINT *,'  END OF INPUT DATA '     !Rev 7/13/90
      CALL ENDPLT
33    IF(PLT.NE.0.d0) CALL DDEND
      bell='16#07'                 !Rev 7/13/90 here thru STOP:
      face='16#02'                 ! to display elapsed time and 
      call systim(tod)           ! time completed.
      call syssec(sec) 
      elsec=(sec-stsec)
      if(sec.lt.stsec)elsec=(sec-stsec+86400.)  !add 24 hr
      elmin=elsec/60
      elhrs=elmin/60
      write (*,886)int(elhrs),int(mod(elmin,60.0d0)),mod(elsec,60.0d0)
                                                                   ! 10 Sep 2008
886   format(/'  END OF INPUT DATA...',22x,'ELAPSED TIME = ',
     1  i2,' HR ',i2,' MIN ',f4.1,' SEC')
      write (*,888)bell,face,face,face,tod,face,face,face  
888   format (//12x,a1,a1,1x,a1,1x,a1,
     1   '   HARPO RUN COMPLETED NORMALLY AT ',A8,2x,a1,1x,a1,1x,a1)
      STOP '  '
C
4     READ(NRYIND,'(A)',END=10001) STMP1
C
c      print *, ' STMP1 = ',STMP1
      READ(STMP1,1100) NW,WWW,LINEX
1100  FORMAT (BZ,I3,E14.7,A)
      IF(NDEVTMP.GT.0.AND.(NFRMAT.NE.-2.OR.NW.LT.0))
     1      WRITE(NDEVTMP,'(A)') STMP1
10001 IF (NW.EQ.0) GO TO 10
      IF(NW.GT.MAXW) GO TO 3400
C
      IF (NW.LE.MAXW .AND. NW.GT.0) GO TO 5
C     TABULAR INPUT DATA
      NWOK=.FALSE.
C 'OPEN' THE TEMP FILE
1150  FORMAT(I3,T18,A)
      NWP=-NW
      IF(NFRMAT.EQ.-2) THEN
C     USE DUMMY ARRAYS TO ABSORB TABULAR DATA
         CALL READW1(AB,NWOK,NW,XMX,XNTBL,XITBL,XFRMTBL,XGP)
      ELSEIF(NWP.EQ.1) THEN
         CALL READW1(AB,NWOK,NW,UMX,UNTBL,UITBL,UFRMTBL,UGP)
      ELSEIF(NWP.EQ.2) THEN
         CALL READW1(AB,NWOK,NW,DUMX,DUNTBL,DUITBL,DUFRMTB,DUGP)
      ELSEIF(NWP.EQ.3) THEN
         CALL READW1(AB,NWOK,NW,CMX,CNTBL,CITBL,CFRMTBL,CGP)
      ELSEIF(NWP.EQ.4) THEN
         CALL READW1(AB,NWOK,NW,DCMX,DCNTBL,DCITBL,DCFRMTB,DCGP)
      ELSEIF(NWP.EQ.5) THEN
         CALL READW1(AB,NWOK,NW,TMX,TNTBL,TITBL,TFRMTBL,TGP)
      ELSEIF(NWP.EQ.6) THEN
         CALL READW1(AB,NWOK,NW,DTMX,DTNTBL,DTITBL,DTFRMTB,DTGP)
      ELSEIF(NWP.EQ.7) THEN
         CALL READW1(AB,NWOK,NW,MMX,MNTBL,MITBL,MFRMTBL,MGP)
      ELSEIF(NWP.EQ.8) THEN
         CALL READW1(AB,NWOK,NW,RMX,RNTBL,RITBL,RFRMTBL,RGP)
      ELSEIF(NWP.EQ.9) THEN
         CALL READW1(AB,NWOK,NW,GMX,GNTBL,GITBL,GFRMTBL,GGP)
      ELSEIF(NWP.EQ.10) THEN
         CALL READW1(AB,NWOK,NW,DGMX,DGNTBL,DGITBL,DGFRMTB,DGGP)
      ELSEIF(NWP.EQ.11) THEN
         CALL READW1(AB,NWOK,NW,SMX,SNTBL,SITBL,SFRMTBL,SGP)
      ELSEIF(NWP.EQ.12) THEN
         CALL READW1(AB,NWOK,NW,DSMX,DSNTBL,DSITBL,DSFRMTB,DSGP)
      ELSEIF(NWP.EQ.13) THEN ! 31 Jan 2006
         CALL READW1(AB,NWOK,NW,AMX,ANTBL,AITBL,AFRMTBL,AGP) ! 31 Jan 2006
      ELSEIF(NWP.EQ.14) THEN ! 31 Jan 2006
         CALL READW1(AB,NWOK,NW,DAMX,DANTBL,DAITBL,DAFRMTB,DAGP) ! 31 Jan 2006
      ELSEIF(NWP.EQ.15) THEN ! 31 Jan 2006
         CALL READW1(AB,NWOK,NW,DMX,DNTBL,DITBL,DFRMTBL,DGP) ! 31 Jan 2006
      ELSEIF(NWP.EQ.16) THEN ! 31 Jan 2006
         CALL READW1(AB,NWOK,NW,DDMX,DDNTBL,DDITBL,DDFRMTB,DDGP) ! 31 Jan 2006
      ELSEIF(NWP.EQ.17) THEN
         CALL READW1(AB,NWOK,NW,VMX,VNTBL,VITBL,VFRMTBL,VGP)
      ELSEIF(NWP.EQ.18) THEN
         CALL READW1(AB,NWOK,NW,DVMX,DVNTBL,DVITBL,DVFRMTB,DVGP)
      ELSEIF(NWP.EQ.19) THEN
         CALL READW1(AB,NWOK,NW,PRMX,PRNTBL,PRITBL,PRFRMTB,PRGP)
      ELSEIF(NWP.EQ.20) THEN
         CALL READW1(AB,NWOK,NW,DPMX,DPNTBL,DPITBL,DPFRMTB,DPGP)
      ELSEIF(NWP.EQ.21) THEN
         CALL READW1(AB,NWOK,NW,MX21,NTBL21,ITBL21,FRMTB21,GP21)  ! 15 Mar 2010
      ELSEIF(NWP.EQ.22) THEN
         CALL READW1(AB,NWOK,NW,MX22,NTBL22,ITBL22,FRMTB22,GP22)  ! 15 Mar 2010
      ELSEIF(NWP.EQ.23) THEN
         CALL READW1(AB,NWOK,NW,MX23,NTBL23,ITBL23,FRMTB23,GP23)  ! 15 Mar 2010
      ELSEIF(NWP.EQ.24) THEN
         CALL READW1(AB,NWOK,NW,MX24,NTBL24,ITBL24,FRMTB24,GP24)  ! 15 Mar 2010
      ELSE
         CALL SREAD1(AB,NWOK,NW)
      ENDIF
      IF(NWOK) GO TO 4
C
3400  WRITE(3,4000) NW,MAXW
4000  FORMAT(15H1THE SUBSCRIPT ,I3,  ' ON THE W-ARRAY INPUT IS OUT OF BO
     1UNDS.  ALLOWABLE VALUES ARE -8 THROUGH,',I4,'.' )
      STOP 4000
C
5     READ(LINEX,70) DEG,KM,NM,FEET,CYCLE,PER,MSKMS
70    FORMAT(7A)
      W(NW)=WWW
C
C     CHECK FOR A 'TERRAIN RELATIVE' HEIGHT SPECIFICATION.
C**   IF SO ADD FLAG VALUE 10**40 TO BE TESTED FOR LATER IN 'RAYTRC'
C     IF SO ADD FLAG VALUE 10**25 TO BE TESTED FOR LATER IN 'RAYTRC' 28 Nov 2003
C
C**   IF(MSKMS.EQ.'T' .AND. WWW.EQ.0.0) W(NW)=1.E-40
c     IF(MSKMS.EQ.'T' .AND. WWW.EQ.0.d0) W(NW)=1.d-10 
      IF(MSKMS.EQ.'T' .AND. WWW.EQ.0.d0) W(NW)=1.d-25 ! 28 Nov 2003
      IF(WWW.EQ.0.d0) GO TO 4
C**   IF(MSKMS.EQ.'T') WWW=WWW*1.E40
c     IF(MSKMS.EQ.'T') WWW=WWW*1.d10
      IF(MSKMS.EQ.'T') WWW=WWW*1.d25                  ! 28 Nov 2003
C
C     CHECK FOR KEYWORD UNITS SPECIFICATION, IF SO PERFORM CONVERSION
C
      uucon=  UCON(KU,LINEX(:10),CONV)
      if (.not. uucon) go to 60
75    IF(CONV.GE.0.d0) THEN
C         W(NW)=WWW*SNGL(CONV)                         ! 5 January 2006
         W(NW)=WWW*CONV
      ELSEIF(W(NW).NE.0.d0) THEN
         W(NW)=-CONV/W(NW)
      ENDIF
      GO TO 4
C
60    IF (DEG.EQ.'1') WWW=WWW*RAD
      IF (KM.EQ.'1') WWW=WWW/EARTHR
      IF (NM.EQ.'1') WWW=WWW*1.852d0
      IF (FEET.EQ.'1') WWW=WWW*3.048006096d-4
      IF(CYCLE.EQ.'1') WWW=WWW*PIT2
      IF(PER.EQ.'1') WWW=PIT2/WWW
      IF(MSKMS.EQ.'1') WWW=WWW*1.d-3
      W(NW)=WWW
      GO TO 4
C
10    continue
cc      print* , ' leaving readw because AB = ', AB
      IF(.NOT.AB) RETURN
      PRINT 1200
1200  FORMAT(' A DATA FORMATTING ERROR PREVENTS CONTINUED EXECUTION')
      STOP 4500
C
C ----------------------------------------------------------------------!-------
      ENTRY SETW(XX)
C ----------------------------------------------------------------------!-------
C     THIS ENTRY IS CALLED ONCE BEFORE THE FIRST RUN SET IS READ
C
C INITIALIZE SOME MATHEMATICAL CONSTANTS
      PI=4.d0*dATAN(1.d0)
      PIT2=2.d0*PI
      PID2=PI/2.d0
      DEGS=180.d0/PI
      RAD=PI/180.d0
      ALN10=dLOG(10.d0)
CC********* INITIALIZE SOME VARIABLES IN THE W ARRAY
      MAXW=1000
      CALL CLEAR(W,MAXW)
      PLON=0.d0
      PLAT=PID2
      EARTHR=6370.d0
      INTYP=3.d0
      MAXERR=1.d-4
      ERATIO=50.d0
      STEP1=1.d0
      STPMAX=100.d0
      STPMIN=1.d-8
      FACTR=5.d-1
      HITLET=.15d0
      MAXSTP=1000.d0
      HOP=1.d0
      HMAX=500.d0
      EXTINC=999.999d0
      PAGLN=66.d0
C
C%      CALL SETUCON
      uucon=setucon(x,y,z)
      setw = 1.0         ! 6 Mar 2006
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY SETOVR (XX)
C ----------------------------------------------------------------------!-------
C     THIS ENTRY IS CALLED AFTER EACH RUN SET IS READ
C     PERFORM OVERRIDE ASSIGNMENTS TO W-ELEMENTS
      CALL OVERRD(MAXSTP,0.d0,1.d3,WCOMTS(23),2,1)
      CALL OVERRD(ERATIO,0.d0,5.d1,WCOMTS(43),2,1)
      CALL OVERRD(FACTR,0.d0,5.d-1,WCOMTS(47),2,1)
	bigval = I1mach(9)                                         ! 26 Aug 2008
      CALL OVERRD(SKIP,0.d0,bigval,WCOMTS(71),2,1)
      CALL OVERRD(HITLET,0.d0,1.5d-1,WCOMTS(75),2,1)
      IF(PAGLN.LT.3.d1) PAGLN=0.d0
      CALL OVERRD(PAGLN,0.d0,6.6d1,WCOMTS(77),2,1)
      CALL OVERRD(PFACTR,0.d0,1.d0,WCOMTS(82),2,1)
      setovr = 1.0         ! 6 Mar 2006
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY PRINTW (XX)
C ----------------------------------------------------------------------!-------
C     THIS ENTRY IS CALLED AFTER EACH RUN SET IS READ
C     TO PRINT VALUES OF THE W-ARRAY IN A FORMAT WHICH SHOWS
C     FULL ACCURACY AND THE EFFECTS OF ANY CONVERSIONS OR OVERRIDES.
C
C     GO TO NEXT PAGE, PUT OUT COPY OF INPUT DATA FOR THIS RUN SET
C     ASSUMING WE HAD A TEMPORARY FILE TO USE
C
      IF(NDEVTMP.LE.0) GO TO 1065
      CALL NEWPAG(NPAG,idint(pagln),PC)
      LNSXPG=LINES+idint(pagln)
      CALL PUTHDR(3,PC,NPAG)
      CALL PUTDVR(3)
      CALL PUTKBK(3,1)
C%      CALL PUTKCT(3,'INPUT DATA FILE FOR RUN SET NUMBER'
C%     1   //NUMSTG(NSET,1,'(I5)'))
      WRITE (3,1055)NSET
1055  FORMAT (47X,'INPUT DATA FILE FOR RUNSET NUMBER ',I3)
      CALL PUTKBK(3,1)
      CALL PUTDVR(3)
      CALL PUTKBK(3,1)
cc      print* , 'ndevtmp = ', ndevtmp
      REWIND NDEVTMP
C
1060  READ(NDEVTMP,'(A)',END=1065) LINEX
      IF(LINES.GE.LNSXPG) THEN
         CALL NEWPAG(NPAG,idint(pagln),PC)
         LNSXPG=LINES+idint(pagln)
         CALL PUTHDR(3,PC,NPAG)
         CALL PUTKBK(3,1)
      ENDIF
      CALL PUTKST(3,' '//LINEX)
      GO TO 1060
1065  CONTINUE
C
C     GO TO NEXT PAGE, PUT OUT W-ELEMENTS SHOWING FULL PRECISION
C     USE A TWO COLUMN FORMAT
C
      CALL NEWPAG(NPAG,idint(pagln),PC)
      CALL PUTHDR(3,PC,NPAG)
      CALL PUTDVR(3)
      CALL PUTKBK(3,1)
      CALL PUTKCT(3,'INITIAL VALUES FOR THE W ARRAY')
      CALL PUTKCT(3,'ONLY NONZERO VALUES PRINTED'
     1  //' -- ALL ANGLES IN RADIANS')
      CALL PUTKBK(3,1)
      CALL PUTDVR(3)
      LINSPP=idint(pagln)-7
C
      WFRMT='(BZ,I3,E14.7,A)'
      IF(NFRMAT.EQ.0) WFRMT='(I4,E24.15,A)'
C
      NWE=0
C     ALLOW MAXIMUM OF 10 PAGES OF OUTPUT(1 IS ENOUGH)
      DO 18 NPGS=1,10
C     ALLOW 3 LINES FOR INSIDE HEADER(SEE BELOW)
         LINSPP=LINSPP-3
C     NEXT ELEMENT TO SCAN IS ONE MORE THAN LAST ONE
         LWE=NWE+1
C
C     INITIALIZE SECOND COLUMN START INDEX
         NXCL=0
C     INITIALIZE TOTAL ELEMENT COUNTER
         NELS=0
C     LOOP TO FIND BREAK POINTS FOR FIRST/SECOND COLUMNS
         DO 14 I=LWE,MAXW
            IF(W(I).EQ.0.d0) GO TO 14
            NWE=I
            NELS=NELS+1
            IF(NELS.EQ.LINSPP) NXCL=I
            IF(NELS.EQ.2*LINSPP) GO TO 15
14       CONTINUE
15    IF(NELS.EQ.0) GO TO 22
C
         IF(NPGS.GT.1) THEN
            CALL NEWPAG(NPAG,idint(pagln),PC)
            CALL PUTHDR(3,PC,NPAG)
            LINSPP=idint(pagln)-1
         ENDIF
C
C     INSERT 'INNER' HEADER
         CALL PUTKBK(3,2)
         STMP1=' N              W(N)'
         CALL PUTKST(3,STMP1(:65)//' '//STMP1)
C
       NX=NXCL-1
       IF(NX.LT.LWE) NX=NWE
       DO 17 NW=LWE,NX
       NCOMT=MIN0(NW,MXCMTS)
       IF(W(NW).NE.0.d0) THEN
         WRITE(STMP1,WFRMT) NW,W(NW),WNOTES(WCOMTS(NCOMT))
         STMP2=' '
         IF(NXCL.EQ.0) GO TO 23
         DO 19 I=NXCL,MAXW
         IF(W(I).NE.0.d0) THEN
           NCOMT=MIN0(I,MXCMTS)
           WRITE(STMP2,WFRMT) I,W(I),WNOTES(WCOMTS(NCOMT))
           NXCL=I+1
           GO TO 23
         ENDIF
19       CONTINUE
23       CALL PUTKST(3,STMP1(:65)//' '//STMP2)
       ENDIF
17     CONTINUE
      IF(NELS.EQ.0) GO TO 22
18    CONTINUE
C
22    CONTINUE
C ----------------------------------------------------------------------!-------
      END FUNCTION   READW 
C ----------------------------------------------------------------------!-------
 
 


C ----------------------------------------------------------------------!-------
C from S1
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE CLEAR(A,N)
C ----------------------------------------------------------------------!-------
C**********************************************************************
C     SET N ELEMENTS OF ARRAY A TO ZERO
 
      implicit double precision (a-h,o-z)
 
      DOUBLE PRECISION A(N)
      IF(N.LE.0) RETURN
      DO    10 I=1,N
10    A(I)=0.d0
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE CLEAR
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      FUNCTION ND2B(INDEC)
C ----------------------------------------------------------------------!-------
C
C     CONVERT A NUMERIC DECIMAL DIGIT STRING TO A BIT STRING
C     WITH EACH BIT SET BY A CORRESPONDING DECIMAL DIGIT.
C
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      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
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      LOGICAL FUNCTION UCON(JC1,U,CONV)
C ----------------------------------------------------------------------!-------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      LOGICAL SETUCON
C     PROVIDES KEYWORD UNITS CONVERSION FOR W-ARRAY INPUT.
C     UNITS SPECIFICATION MUST BE IN THE FORM 'DV UN' WHERE DV IS THE
C     TYPE OF VARIABLE(SUCH AS AN FOR ANGLE) AND UN IS THE UNITS CHOICE
C     SUCH AS RD FOR RADIANS). RETURN VALUE IS THE CONVERSION FACTOR.
C     COMMON DECK "CUCON" INSERTED HERE
 
      COMMON/UCONV/CNVV(5,4)
      CHARACTER PCV*3,CNVC*2
      COMMON/UCONC/PCV(5),CNVC(5,4)
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 "WWR" INSERTED HERE
      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))
      CHARACTER U*(*),PU*3,PV*2
C
      DATA NPV,NPU/5,4/
      UCON=.FALSE.
      CONV=1.0d0
C

      READ(U,100) PU,PV
100   FORMAT(A,A)
C
      DO 2010 J=1,NPV
2010     IF(PU.EQ.PCV(J)) GO TO 2015
      JC1=4
      RETURN
2015  JC1=J
      DO 2020 K=1,NPU
2020     IF(PV.EQ.CNVC(J,K)) GO TO 2025
      JC1=4
      RETURN
2025  UCON=.TRUE.
      CONV=CNVV(JC1,K)
      IF(CONV.LT.0.d0.AND.CONV.GT.-2.0D0) CONV=-CONV/EARTHR    ! 19, 20 Mar 2009
      RETURN
C
C     INITIAL CONVERSION CONSTANTS
C ----------------------------------------------------------------------!-------
      ENTRY SETUCON(x,y,z)
C ----------------------------------------------------------------------!-------
 
C*    A bug??? A Newhall
c*      ENTRY SETUCON
 
      SETUCON=.TRUE.
C
c      print *,'RAD =',RAD
      CNVV(1,2) = RAD
C     USE NEGATIVE TAG FOR DYNAMIC USE OF EARTHR                   ! 19 Mar 2009
      CNVV(1,3) = -1.d0
      CNVV(1,4) = -1.852d0                                         ! 19 Mar 2009
      CNVV(2,2) = 1.0d-3
      CNVV(2,3) = 3.048006096d-4
      CNVV(2,4) = 1.852d0
      CNVV(3,2) = PIT2
      CNVV(3,3) = PIT2*1.0d3
C     USE NEGATIVE TAG FOR ENTERING FREQUENCY AS A PERIOD IN SEC   ! 20 Mar 2009
      CNVV(3,4) = -PIT2
C      CNVV(4,1) = aln10/10.d0
      CNVV(4,1) = aln10/20.d0 ! 23 Feb 2007 
C ----------------------------------------------------------------------!-------
      END FUNCTION UCON
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
C from S2
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      FUNCTION GET1(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(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
C ----------------------------------------------------------------------!-------


 
 
C ----------------------------------------------------------------------!-------
      FUNCTION GET(Z)
C ----------------------------------------------------------------------!-------
 
C     'GET' AND ENTRY 'DOT' PROVIDE A CONTROL METHOD FOR AVOIDING
C     REDUNDANT CALLS TO THE TERRAIN AND RECEIVER MODELS. THE VALUES
C     RETURNED ARE THE FUNCTION VALUES FOR 'F' OR 'G' OR THEIR DERIVATI
C     (VIA 'DOT' ENTRY). UNNECESSARY CALLS ARE ELIMINATED THROUGH USE
C     OF 'TIME OF CALL' VARIABLES WHICH ARE COMPARED WITH THE CURRENT
C     LAST CALL TIME MAINTAINED BY THE 'RKAM' PROGRAM. WHEN VALUES ARE
C     NOT CURRENT THEY ARE UPDATED BY CALLS TO THE APPROPRIATE ROUTINES
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      double precision d1mach                                     ! 10 Nov 2009
      DOUBLE PRECISION Z(12),PF(10),PG(10),pu(10)
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
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)),
     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
      DOUBLE PRECISION PI(10),PJ(10), PK(10) ! 10 Jun 2003
c%      COMMON/TMP1/NC,NSZR    ,NSZRR    ,NSZRTH
c%     1 ,NSZRPH    ,NSZTH    ,NSZPH    ,NSZTHTH    ,NSZPHPH    ,NSZTHPH
c%     2 ,NVELECT    ,NWIME
C
      equivalence (pf,pfr),(pg,pgr),(pu,pur)
      DATA PI,PJ,pk/1.d0,9*0.d0,1.d0,9*0.d0,1.d0,9*0.d0/
c%    DATA NC,NSZR,NSZRR,NSZRTH,NSZRPH,NSZTH,NSZPH,NSZTHTH,NSZPHPH
c%   1   ,NSZTHPH,NVELECT,NWIME/1,2,3,4,5,6,7,8,9,10,11,12/
c%    DATA FWIME,GWIME,FVELECT,GVELECT/2*-1.d0,'RECEIVER','TERRAIN' /
c%    DATA UWIME,UVELECT/-1.d0,'SURFACE' /
c%    DATA ZHRO/0.d0/
C
      IENTRY=1
      GO TO 5
 
C ----------------------------------------------------------------------!-------
      ENTRY SETGET (x)
C ----------------------------------------------------------------------!-------
      CALL RMOVE(PF,PI,10)
      CALL RMOVE(PG,PJ,10)
      call rmove(pu,pk,10)
c%    CALL IMOVE(NZ,NC,12)
      rktime=0.d0
! 25 Aug 2011      fselect='RECEIVER'
! 25 Aug 2011      gselect='TERRAIN'
! 25 Aug 2011      uselect='SURFACE'
      fselect=8HRECEIVER                                           ! 25 Aug 2011
      gselect=7HTERRAIN                                            ! 25 Aug 2011
      uselect=7HSURFACE                                            ! 25 Aug 2011
      ftime=-1.d0
      utime=-1.d0
      gtime=-1.d0
      zero=0.d0
      nz=1
      npzr=2
      npzrr=3
      npzrth=4
      npzrph=5
      npzth=6
      npzph=7
      npzthth=8
      npzphph=9
      npzthph=10
      nselect=11
      ntime=12
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY DOT(Z)
C ----------------------------------------------------------------------!-------
      IENTRY=2
5     CONTINUE
c diagnostic      print *,ftime,gtime,utime
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
      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
      GET=Z(NZ)
      RETURN
C
20    GET=Z(NPZR)*DRDT(1)+Z(NPZTH)*DRDT(2)+Z(NPZPH)*DRDT(3)
      dot=get
      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION GET
C ----------------------------------------------------------------------!-------
 
 
 
 
C ----------------------------------------------------------------------!-------
      FUNCTION ITEST(I)
C ----------------------------------------------------------------------!-------
C     USED TO PASS INTEGER VALUES THROUGH FOR VARIABLES TYPED REAL
C     (AS IN MIXED REAL/INTEGER ARRAYS)
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      ITEST=I
C ----------------------------------------------------------------------!-------
      END FUNCTION ITEST
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      FUNCTION ITOC(N)
C ----------------------------------------------------------------------!-------
C     RETURN 7 CHARACTER STRING REPRESENTATION OF INTEGER N
C     IF NUMBER IS TOO LARGE OR SMALL USE FLOATING POINT FORMAT
      CHARACTER ITOC*7
      IF(N.LT.-9999.OR.N.GT.99999) GO TO 100
      ITOC=' '
      WRITE(ITOC,'(I7)',ERR=100) N
      RETURN
100   WRITE(ITOC,'(2PG7.0)') FLOAT(N)
      return
C ----------------------------------------------------------------------!-------
      END FUNCTION ITOC
C ----------------------------------------------------------------------!-------
 
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE SETTRC
C ----------------------------------------------------------------------!-------
C     INITIALIZE /TRAC/ LABELED COMMON

C     COMMON DECK "TRAC" INSERTED HERE
 
      implicit double precision (a-h,o-z)
 
      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
C ----------------------------------------------------------------------!-------
 
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      LOGICAL FUNCTION WCHANGE(W1,W2)
C ----------------------------------------------------------------------!-------
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      DOUBLE PRECISION W1(400),W2(400)
C
C WCHANGE PROVIDES A TEST AGAINST TWO W-ARRAYS
C     IF BOTH WOULD PRODUCE THE SAME SET OF RAYS BY RAYTRC THEN RESULT
C     IS <FALSE>.
C
      INTEGER NDX(8)
C
      DATA NGRPS,NDX/8, 1,17, 21,26, 41,47, 100,399/
C
      WCHANGE=.FALSE.
      DO 20 N=1,NGRPS,2
         N1=NDX(N)
         N2=NDX(N+1)
      DO 20 I=N1,N2
          WCHANGE=W1(I).NE.W2(I)
          IF(WCHANGE) RETURN
20    CONTINUE
C ----------------------------------------------------------------------!-------
      END FUNCTION WCHANGE
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      FUNCTION RENORM(VECTOR,NNORM,NCOMPS)
C ----------------------------------------------------------------------!-------
C     NORMALIZES 'NCOMPS' COMPONENT VECTOR 'VECTOR' TO MAGNITUDE
C     'NNORM' AND RETURNS SQUARE ROOT OF FACTOR NEEDED.
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
      DOUBLE PRECISION NNORM,VECTOR(*)
C
      RENORM=0.d0
      IF(NNORM.LE.0.d0) RETURN
C
      RENORM=0.d0
      DO 10 I=1,NCOMPS
10    RENORM=RENORM+VECTOR(I)*VECTOR(I)
      IF(RENORM.EQ.0.d0) RETURN
C
      RENORM=DSQRT(NNORM/RENORM)
      DO 20 I=1,NCOMPS
20    VECTOR(I)=VECTOR(I)*RENORM
C
      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION RENORM
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE SET2(A,V,N)
C ----------------------------------------------------------------------!-------
C     SETS N COMPONENTS OF VECTOR A TO SINGLE VALUE V
      DOUBLE PRECISION A(N),v
      INTEGER NA(N)
      CHARACTER NV
      DO 100 I=1,N
100   A(I)=V
      RETURN
C ----------------------------------------------------------------------!-------
	end SUBROUTINE SET2 ! 10 Jun 2003, added this and rewrote below
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE ISET2(NA,NV,N)
C	ENTRY ISET2(NA,NV,N)
C ----------------------------------------------------------------------!-------
C     SPECIAL ENTRY FOR INTEGER ARRAYS FOR MACHINES HAVING DIFFERENT
C     WORD SIZES FOR INTEGERS THAN REALS.
	INTEGER N,NA(N)
      CHARACTER NV
      DO 200 I=1,N
C200	NA(I) = NV
200   NA(I)=IACHAR(NV)      
C200	continue
	RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE ISET2
C ----------------------------------------------------------------------!-------
 
 


C ----------------------------------------------------------------------!-------
C from S3
C ----------------------------------------------------------------------!-------
      SUBROUTINE PUTDES(NUNIT,DES,MOD,ids)
C ----------------------------------------------------------------------!-------
C     TITLING AID, OUTPUT MODULE DESCRIPTION AS HOLLERITH NAME
C     AND NUMERIC IDENTIFIER
      implicit double precision (a-h,o-z)
      CHARACTER DES*(*),TITLE*30,SMODL*20,NUMSTG*80,strg5*200
      DOUBLE PRECISION MOD(2)
      TITLE=DES
      WRITE(SMODL,'(A8,F10.2)')MOD
      CALL PUTKST(NUNIT, ' '//TITLE//SMODL//NUMSTG(ids,20,'(20a4)'))
      return
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PUTDES
C ----------------------------------------------------------------------!-------
 
 
 
 
C ----------------------------------------------------------------------!-------
      FUNCTION NUMSTG(V,N,FRM)
C ----------------------------------------------------------------------!-------
C     CONVERTS AN integer v(n) to a string numstg*80
C     USING FORMAT frm PROVIDED IN CALL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER NUMSTG*80,FRM*(*)
      INTEGER V(N)
C
      NUMSTG=' '
      WRITE(NUMSTG,FRM) V
      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION NUMSTG
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE SFILL(STG,LN,C)
C ----------------------------------------------------------------------!-------
C     FILLS A STRING WITH N SPECIFIED CHARACTERS
      implicit double precision (a-h,o-z)
      CHARACTER STG*(*),C*1
C
      DO 10 I=1,LN
10    STG(I:I)=C
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE SFILL
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      INTEGER FUNCTION STRIM(C)
C ----------------------------------------------------------------------!-------
C     DETERMINES POSITION OF LAST NONBLANK CHARACTER OF A STRING
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER*(*) C
C
      L=LEN(C)
      DO 100 I=L,1,-1
100   IF(C(I:I) .NE. ' ') GO TO 200
      I=0
200   STRIM=I
      return
C ----------------------------------------------------------------------!-------
      END FUNCTION STRIM
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      FUNCTION RERR(V,ERT,ELAB,NKV,PREF)
C ----------------------------------------------------------------------!-------
C     RETURNS RELATIVE ERROR OF VARIABLE RKVAR(NKV) IF PREF<>0 AND
C     THE ERROR IS GREATER THAN PREVIOUS ERROR 'V'.
C     COMMON DECK "RKAM" INSERTED HERE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      character*(*) elab
      DOUBLE PRECISION KR,KTH,KPH
      COMMON//R,TH,PH,KR,KTH,KPH,RKVARS(14),TPULSE,CSTEP,DRDT(20)
C
      RERR=V
      IF(NKV.LE.0) RETURN
      IF(PREF.EQ.0.d0) RETURN
      V1=RKVARS(NKV)/PREF-1.d0
      IF(dABS(V1).LT.dABS(V)) RETURN
C
      ERT=IACHAR(ELAB)
      RERR=V1
C ----------------------------------------------------------------------!-------
      END FUNCTION RERR
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE RERROR(ROUTIN,STR,VAL)
C ----------------------------------------------------------------------!-------
C     REPORTS ERROR CONDITIONS AND STOPS PROGRAM.
C
C%      implicit double precision (a-h,o-z)
	character routin*(*),str*(*)
      PRINT 10,ROUTIN,STR,VAL
10    FORMAT(' ERROR CONDITION IN RAYTRACE ROUTINE <'
     1  ,A6,'> DUE TO ',A10,/,10x,' VALUE=  ',F8.2)
      STOP
C ----------------------------------------------------------------------!-------
      END SUBROUTINE RERROR
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE STOPIT(A)
C ----------------------------------------------------------------------!-------
C     PRINTS CONDITION AND STOPS PROGRAM
C     AFTER CALLING THE SYSTEM POST MORTEM DUMP.
      implicit double precision (a-h,o-z)
      CHARACTER A*(*)
C
      PRINT 100, A
      CALL MORTEM
100   FORMAT('*** STOPIT WITH CONDITION <',A,'>')
      STOP
C ----------------------------------------------------------------------!-------
      END SUBROUTINE STOPIT
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE PUTKST(NUNIT,STRG)
C ----------------------------------------------------------------------!-------
C     WRITE LINE OF OUTPUT TO PRINTER UNIT ADDING TO LINE COUNT
      implicit double precision (a-h,o-z)
      CHARACTER STRG*(*)
      INTEGER STRIM
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
C     PUT OUT A STRING WITH LINE COUNT INCREMENT
      LINES=LINES+1
      LN=MAX0(1,MIN0(STRIM(STRG),132))
      WRITE(NUNIT,'(A)') STRG(1:LN)
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PUTKST
C ----------------------------------------------------------------------!-------
C
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE PUTKCT(NUNIT,STRG)
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)
C     PUT OUT A CENTERED LINE AND COUNT
      CHARACTER STRG*(*)
      CHARACTER BLANKS*100,STMP*80
      INTEGER STRIM
      CHARACTER pcc*(*),PC*1,typ*20
C
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
      DATA BLANKS/' '/
C
      NTRM=MAX0(1,STRIM(STRG))
      NBLKS=66-(NTRM+1)/2
      LINES=LINES+1
      LN=MIN0(NTRM,132-NBLKS)
      STMP=STRG
      WRITE(NUNIT,'(A)') BLANKS(:NBLKS)//STMP(:LN)
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY PUTKBK(NUNIT,NBKS)
C ----------------------------------------------------------------------!-------
C     PUT NBLKS BLANK LINES TO UNIT NUNIT
C
      LINES=LINES+1
      DO 100 I=1,NBKS
100      WRITE(NUNIT,'(1X)')
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY NEWPAG(NXP,LIP,PCC)
C ----------------------------------------------------------------------!-------
C............NXP= page number (input)
C............LIP= lines per page (input)
C............PC= new page control character (output)
C     THIS ENTRY COMPUTES NEXT PAGE NUMBER AND LINE COUNT TO END
C     OF CURRENT PAGE
C     RETAIN LINES PER PAGE FOR USE BY TSTPAG ENTRY
      LINSPP=LIP
C     THE LOGIC HERE ASSUMES THAT A FORM FEED IS ALWAYS NEEDED.
C     A MORE SOPHISTICATED VERSION WOULD DETERMINE IF A NEW PAGE
C     IS NEEDED AND RETURN PC=' ' IF NOT. THE CALLING PROGRAM COULD
C     USE OR IGNORE.
C%    PC='1'
      call pcodes (pcc,typ)
C     COMPUTE CURRENT PAGE NUMBER
      NPAG=(LINES+LINSPP-1)/LINSPP
C     COMPUTE LINE COUNT AT END OF CURRENT PAGE
      LINES=NPAG*LINSPP
      NXP=NPAG+1
C     COMPUTE LINE COUNT AT END OF NEXT PAGE
      LNSXPG=LINES+LINSPP-1
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY TSTPAG(NUNIT)
C ----------------------------------------------------------------------!-------
C     OUTPUT HEADER IF HAVE REACHED LAST LINE OF PAGE
      IF(LINES.GE.LNSXPG) THEN
C     COMPUTE CURENT PAGE NUMBER
         NPAG=(LINES+LINSPP-1)/LINSPP
C     AND ITS LAST LINE
         LINES=NPAG*LINSPP
         LNSXPG=LINES+LINSPP-1
         NPAG=NPAG+1
         call pcodes (pc,typ)
         CALL PUTHDR(NUNIT,pc,NPAG)
         WRITE(NUNIT,'(1X)')
         LINES=LINES+1
      ENDIF
C ----------------------------------------------------------------------!-------
      END SUBROUTINE PUTKCT
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE OPNREP(IUN,FNAME)
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)
C     THIS OPEN OPERATION ALLOWS FOR AN EXISTING VERSION OF
C     A FILE. IF THE FILE EXISTS IT IS OVERWRITTEN.
C     THE LOGIC USED IS HOPEFULLY MACHINE INDEPENDENT.
      CHARACTER*(*) FNAME
      CALL ZAPFIL(IUN,FNAME)
      OPEN(IUN,FILE=FNAME,STATUS='NEW')
c      print* , ' iun= ', iun
      REWIND IUN
      RETURN
C ----------------------------------------------------------------------!-------
      ENTRY OPNURP(IUN,FNAME)
C ----------------------------------------------------------------------!-------
      CALL ZAPFIL(IUN,FNAME)
      OPEN(IUN,FILE=FNAME,STATUS='NEW',FORM='UNFORMATTED')
c      print* , ' iun= ', iun
      REWIND IUN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE OPNREP
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE ZAPFIL(IUN,FNAME)
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)
C     TEST FOR EXISTENCE OF FILE AND DELETE
      CHARACTER FNAME*(*)
      LOGICAL EX
      INQUIRE(FILE=FNAME,EXIST=EX)
      IF(EX) THEN
         OPEN(UNIT=IUN,FILE=FNAME,STATUS='OLD')
         CLOSE(UNIT=IUN,STATUS='DELETE')
      ENDIF
C ----------------------------------------------------------------------!-------
      END SUBROUTINE ZAPFIL
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE OVERRD(VAR,TST,DEFT,NFLG,NFVLEQ,NFVLNE)
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)
C     OVERRIDE SUPPORT ROUTINE
C     TEST 'VAR' AGAINST 'TST', IF EQUAL SET TO DEFAULT 'DEFT'
C     AND ALSO SET INTEGER FLAG 'NFLG' TO VALUE 'NFVLEQ' ELSE 'NFVLNE'
C
      IF(VAR.EQ.TST) THEN
         VAR=DEFT
         NFLG=NFVLEQ
      ELSE
         NFLG=NFVLNE
      ENDIF
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE OVERRD
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE SFILTR(C,S,KSET)
C ----------------------------------------------------------------------!-------
C     FILTERS EXTRANEOUS CHARACTERS FROM PLOT LABELS
      implicit double precision (a-h,o-z)
      CHARACTER*(*) C,S,KSET
C
      LN=LEN(C)
      J=0
      DO 10 I=1,LN
        IF(INDEX(KSET,C(I:I)).EQ.0) THEN
           J=J+1
           S(J:J)=C(I:I)
        ENDIF
10    CONTINUE
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE SFILTR
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      FUNCTION ALCOSH(X)
C ----------------------------------------------------------------------!-------
C     COMPUTE LOG(COSH(X)) AND USE LARGE ARGUMENT APPROXIMATION
C     WHEN POSSIBLE.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DATA ALOG2/.6931471805599453d0/
C
      IF(dABS(X).GT.50.d0) GO TO 10
      EX=dEXP(X)
      ALCOSH=DLOG((EX+1.D0/EX)*.5D0)
      RETURN
10    ALCOSH=DABS(X)-alog2
C ----------------------------------------------------------------------!-------
      END FUNCTION ALCOSH
C ----------------------------------------------------------------------!-------
 
 
 
 
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE GAUSEL (C,NRD,NRR,NCC,NSF)
C ----------------------------------------------------------------------!-------
C     MATRIX INVERSION BY METHOD OF GAUSSIAN PIVOT
C     SOLVES NCC-NRR SETS OF NRR EQUATIONS
C********* SAME AS SUBROUTINE GAUSSEL WRITTEN BY L. DAVID LEWIS ********
      implicit double precision (a-h,o-z)
      DIMENSION C(NRD,NCC),L(128,2)
C          BITS = 2.**-18
      DATA BITS/3.8146972656d-6/
      NR=NRR
      NC=NCC
      IF(NC.LT.NR.OR.NR.GT.128.OR.NR.LE.0) STOP 5000
C
C          INITIALIZE.
      NSF=0
      NRM=NR-1
      NRP=NR+1
      D=1.
      LSD=1
      DO 1 KR=1,NR
      L(KR,1)=KR
1     L(KR,2)=0
      IF(NR.EQ.1) GO TO 42
C
C          ELIMINATION PHASE.
      DO 41 KP=1,NRM
      KPP=KP+1
      PM=0.d0
      MPN=0
C
C          SEARCH COLUMN KP FROM DIAGONAL DOWN FOR MAX PIVOT.
      DO 2 KR=KP,NR
      LKR=L(KR,1)
      PT=dABS(C(LKR,KP))
      IF(PT.LE.PM) GO TO 2
      PM=PT
      MPN=KR
      LMP=LKR
2     CONTINUE
C
C          IF MAX PIVOT IS ZERO, MATRIX IS SINGULAR.
      IF(MPN.EQ.0) GO TO 9
      NSF=NSF+1
      IF(MPN.EQ.KP) GO TO 3
C
C          NEW ROW NUMBER KP HAS MAX PIVOT.
      LSD=-LSD
      L(KP,2)=L(KP,1)
      L(MPN,1)=L(KP,1)
      L(KP,1)=LMP
C
C          ROW OPERATIONS TO ZERO COLUMN KP BELOW DIAGONAL.
3     MKP=L(KP,1)
      P=C(MKP,KP)
      D=D*P
      DO 41 KR=KPP,NR
      MKR=L(KR,1)
      Q=C(MKR,KP)/P
      IF(Q.EQ.0.d0) GO TO 41
C
C          SUBTRACT Q * PIVOT ROW FROM ROW KR.
      DO 4 LC=KPP,NC
      R=Q*C(MKP,LC)
      C(MKR,LC)=C(MKR,LC)-R
4     IF(dABS(C(MKR,LC)).LT.dABS(R)*BITS) C(MKR,LC)=0.d0
41    CONTINUE
C
C          LOWER RIGHT HAND CORNER.
42    LNR=L(NR,1)
      P=C(LNR,NR)
      IF(P.EQ.0.d0) GO TO 9
      NSF=NSF+1
      D=D*P*LSD
      IF(NR.EQ.NC) GO TO 8
C
C          BACK SOLUTION PHASE.
      DO 61 MC=NRP,NC
      C(LNR,MC)=C(LNR,MC)/P
      IF(NR.EQ.1) GO TO 61
      DO 6 LL=1,NRM
      KR=NR-LL
      MR=L(KR,1)
      KRP=KR+1
      DO 5 MS=KRP,NR
      LMS=L(MS,1)
      R=C(MR,MS)*C(LMS,MC)
      C(MR,MC)=C(MR,MC)-R
5     IF(dABS(C(MR,MC)).LT.dABS(R)*BITS) C(MR,MC)=0.d0
6     C(MR,MC)=C(MR,MC)/C(MR,KR)
61    CONTINUE
C
C          SHUFFLE SOLUTION ROWS BACK TO NATURAL ORDER.
      DO 71 LL=1,NRM
      KR=NR-LL
      MKR=L(KR,2)
      IF(MKR.EQ.0) GO TO 71
      MKP=L(KR,1)
      DO 7 LC=NRP,NC
      Q=C(MKR,LC)
      C(MKR,LC)=C(MKP,LC)
7     C(MKP,LC)=Q
71    CONTINUE
C
C          NORMAL AND SINGULAR RETURNS.  GOOD SOLUTION COULD HAVE D=0.
8     C(1,1)=D
      GO TO 91
9     C(1,1)=0.d0
   91 RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE GAUSEL 
C ----------------------------------------------------------------------!-------
 



C ----------------------------------------------------------------------!-------
      SUBROUTINE FTANH(H,T,PTR,N,Z0,Z,T0,DL,C,CAV,DLC,ALC,DC)
C ----------------------------------------------------------------------!-------
c     28 April 2004
c     FTANH is a general purpose hyperbolic tangent subroutine.
C     TO MAKE A PROFILE REPRESENTED BY A SEQUENCE OF LINEAR SEGMENTS
C     SMOOTHLY JOINED BY HYPERBOLIC FUNCTIONS. PARAMETERS ARE INPUT
C     AS TABULAR DATA WITH SLOPES COMPUTED FROM TABULAR DATA.
C     REFERENCE VALUE T0 IS READ FROM TABULAR DATA.
c     Input: H = independent variable, often height. 
c     Output: T = value of dependent variable, 
c             PTR = derivative of T with respect to H.
c
c     Other variables are input, and give the profile.  
c     They are defined at ENTRY ITANH.
c
C     Compile with options /c (compile only) and /4Yb (error tracing)
      implicit double precision (a-h,o-z)
C.......The changes from the Cyber version indicated in lower case 
C.......below enlarge the data arrays to accommodate 100 profile points.
      DOUBLE PRECISION c(100), tm(99), z(99), dl(99), DLC(100)
     1 , ALC(100), DC(100)

      SUM = 0.d0

C     LOOP TO SUM OVER ALL COEFFICIENTS
C     USE SPECIAL FUNCTION 'ALCOSH' WHICH ALLOWS FOR LARGE ARGUMENTS.
      DO 1 I = 1, N
   1  SUM = SUM + DLC(I) *(ALCOSH((H - Z(I)) / DL(I)) - ALC(I) )
      T = T0 + SUM + CAV * (H - Z0) 
      SUM = 0.d0
      DO 2 I = 1, N
   2  SUM = SUM + DC(I) * (1.d0 + dTANH ((H - Z(I)) /DL(I))) ! 17 Sep 2004
      PTR = C(1) + SUM

      RETURN

      ENTRY ITANH(N,Z0,Z,T0,TM,DL,C,CAV,DLC,ALC,DC)
C
C     CONVERT 'T' ARRAY INPUT(OVERLAYS 'C' ARRAY) TO 'C' ARRAY
c     Input: 
c            N+2 = number of profile points 
c            Z0 = the independent variable at the lowest profile point
c            Z = array of independent variables in profile
c            T0 = the dependent variable at the lowest profile point
c            TM = array of dependent variable in profile
c            DL = array of transition lengths in profile
c     Output:
c            C = array of slopes for profile
C            CAV = the average of the first and last values in the C array
C            DLC, ALC, and DC are arrays that are based on the profiles and
C                 and are calculated once to avoid recalculating each time.
c     Note: we allow the TM and C arrays to overlap to save storage space, 
c           so that TM(1) is the same as C(2), but we make no use of that here.
c           The TM array is no longer needed after the C array is calculated.
c
      TIM1=T0
      ZIM1=Z0
      NP1=N+1
C
      DO 10 I=1,NP1
         TI=TM(I)
         ZI=Z(I)
         C(I)=(TI-TIM1)/(ZI-ZIM1)
         TIM1=TI
10       ZIM1=ZI
C
      CAV = 0.5d0 * (C(1) + C(N + 1)) 
C
      DO 20 I = 1,N
      DLC(I) = DL(I) * (C(I + 1) - C(I)) / 2.d0 
      ALC(I) = ALCOSH((Z(I)-Z0) / DL(I))
      DC(I) = (C(I + 1) - C(I)) / 2.d0 
20    CONTINUE

      RETURN

C ----------------------------------------------------------------------!-------
      END SUBROUTINE FTANH
C ----------------------------------------------------------------------!-------



C ----------------------------------------------------------------------!-------
      SUBROUTINE FTANH2(H,T,PTR,PTRR,N,Z0,Z,T0,DL,C,CAV,DLC,ALC,DC)
C ----------------------------------------------------------------------!-------
c     3 January 2006
c     FTANH2 is a general purpose hyperbolic tangent subroutine.
C     TO MAKE A PROFILE REPRESENTED BY A SEQUENCE OF LINEAR SEGMENTS
C     SMOOTHLY JOINED BY HYPERBOLIC FUNCTIONS. PARAMETERS ARE INPUT
C     AS TABULAR DATA WITH SLOPES COMPUTED FROM TABULAR DATA.
C     REFERENCE VALUE T0 IS READ FROM TABULAR DATA.
C     Input: H = independent variable, often height. 
C     Output: T = value of dependent variable, 
C             PTR = derivative of T with respect to H.
C             PTRR = second derivative of T with respect to H.
c
c     Other variables are input, and give the profile.  
c     They are defined at ENTRY ITANH2.
c
C     Compile with options /c (compile only) and /4Yb (error tracing)
      implicit double precision (a-h,o-z)
C.......The changes from the Cyber version indicated in lower case 
C.......below enlarge the data arrays to accommodate 100 profile points.
      REAL(KIND=8) c(100), tm(99), z(99), dl(99), DLC(100),
     1 ALC(100), DC(100)

C     LOOP TO SUM OVER ALL COEFFICIENTS
C     USE SPECIAL FUNCTION 'ALCOSH' WHICH ALLOWS FOR LARGE ARGUMENTS.
      SUM = 0.d0
      DO 1 I = 1, N
1     SUM = SUM + DLC(I) * (ALCOSH((H - Z(I)) / DL(I)) - ALC(I) )
      T = T0 + SUM + CAV * (H - Z0) 

      SUM = 0.d0
      DO 2 I = 1, N
2     SUM = SUM + DC(I) * (1.d0 + dTANH ((H - Z(I))/DL(I)))
      PTR = C(1) + SUM

      SUM = 0.d0
      DO 3 I = 1, N
3     SUM = SUM + (DC(I) * (1.d0 - (dTANH ((H - Z(I))/DL(I)))**2))
     & * (1/DL(I))
      PTRR = SUM

      RETURN

      ENTRY ITANH2(N,Z0,Z,T0,TM,DL,C,CAV,DLC,ALC,DC)
C
C     CONVERT 'T' ARRAY INPUT(OVERLAYS 'C' ARRAY) TO 'C' ARRAY
c     Input: 
c            N+2 = number of profile points 
c            Z0 = the independent variable at the lowest profile point
c            Z = array of independent variables in profile
c            T0 = the dependent variable at the lowest profile point
c            TM = array of dependent variable in profile
c            DL = array of transition lengths in profile
c     Output:
c            C = array of slopes for profile
C            CAV = the average of the first and last values in the C array
C            DLC, ALC, and DC are arrays that are based on the profiles and
C                 and are calculated once to avoid recalculating each time.
c     Note: we allow the TM and C arrays to overlap to save storage space, 
c           so that TM(1) is the same as C(2), but we make no use of that here.
c           The TM array is no longer needed after the C array is calculated.
c
      TIM1=T0
      ZIM1=Z0
      NP1=N+1
C
      DO 10 I=1,NP1
         TI=TM(I)
         ZI=Z(I)
         C(I)=(TI-TIM1)/(ZI-ZIM1)
         TIM1=TI
10       ZIM1=ZI
C
      CAV = 0.5d0 * (C(1) + C(N + 1)) 
C
      DO 20 I = 1,N
      DLC(I) = DL(I) * (C(I + 1) - C(I)) / 2.d0 
      ALC(I) = ALCOSH((Z(I)-Z0) / DL(I))
      DC(I) = (C(I + 1) - C(I)) / 2.d0 
20    CONTINUE

      RETURN

C ----------------------------------------------------------------------!-------
      END SUBROUTINE FTANH2
C ----------------------------------------------------------------------!-------




C ----------------------------------------------------------------------!-------
C later additions
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE FTANH3(H,T,PTR,PTZ,N,Z0,Z,T0,DL,C,CAV,DLC,ALC,DC
     & ,TC,PCZ)
C ----------------------------------------------------------------------!-------
c     23 Dec 2010
c     FTANH3 is a general purpose hyperbolic tangent subroutine.
C     TO MAKE A PROFILE REPRESENTED BY A SEQUENCE OF LINEAR SEGMENTS
C     SMOOTHLY JOINED BY HYPERBOLIC FUNCTIONS. PARAMETERS ARE INPUT
C     AS TABULAR DATA WITH SLOPES COMPUTED FROM TABULAR DATA.
C     REFERENCE VALUE T0 IS READ FROM TABULAR DATA.
!     FTANH3 is the same as FTANH, except that it also furnishes derivatives of
!     T with respect to z0 and the elements of z
c     Input: H = independent variable, often height. 
c     Output: T = value of dependent variable, 
c             PTR = derivative of T with respect to H.
!             pTz = partial derivatives of T with respect to z(i)
c
c     Other variables are input, and give the profile.  
c     They are defined at ENTRY ITANH3.
c
      implicit double precision (a-h,o-z)
C.......The changes from the Cyber version indicated in lower case 
C.......below enlarge the data arrays to accommodate 100 profile points.
      DOUBLE PRECISION c(100), tm(99), z(99), dl(99), DLC(100)
     1 , ALC(100), DC(100), TC(100), TA(100), AL(100)
     2  ,pTz(0:101),pcz(0:101)

      SUM = 0.d0

C     LOOP TO SUM OVER ALL COEFFICIENTS
C     USE SPECIAL FUNCTION 'ALCOSH' WHICH ALLOWS FOR LARGE ARGUMENTS.
      DO 1 I = 1, N
      AL(I) = ALCOSH((H - Z(I)) / DL(I)) 
      SUM = SUM + DLC(I) *(AL(I) - ALC(I) )
   1  CONTINUE
      T = T0 + SUM + CAV * (H - Z0) 
      SUM = 0.d0
      DO 2 I = 1, N
      TA(I) = dTANH ((H - Z(I)) /DL(I))
      SUM = SUM + DC(I) * (1.d0 + TA(I))                            
   2  CONTINUE
      PTR = C(1) + SUM

      IF(n.LT.1) RETURN
      PTZ(0)=-CAV-(h-z0)*PCZ(1)/2.D0 +DL(1)*PCZ(1)*(AL(1)-ALC(1))/2.D0
     & +DC(1)*TC(1)
      i=n+1
      PTZ(n+1)=
     &           (h-z0)*PCZ(n+1)/2.d0+DL(n)*PCZ(n+1)*(AL(n)-ALC(n))/2.d0
      IF (n.eq.1) PTZ(1)=
     & ((h-z0)*(PCZ(1)-PCZ(2))
     & -DL(1)*(PCZ(2)+PCZ(1))*(AL(1)-ALC(1)))/2.d0
      IF (n.gt.1) THEN
         PTZ(1)=
     &           (h-Z0)*PCZ(1)/2.d0 
     &          -DL(1)*(PCZ(2)+PCZ(1))*(AL(1)-ALC(1))/2.d0
     &          +DL(2)*PCZ(2)*(AL(2)-ALC(2))/2.d0
     &          -DC(1)*(TA(1)+TC(1))
         PTZ(n)=
     &          -(h-Z0)*PCZ(n+1)/2.d0 
     &          -DL(n)*(PCZ(n+1)+PCZ(n))*(AL(n)-ALC(n))/2.d0
     &          +DL(n-1)*PCZ(n)*(AL(n-1)-ALC(n-1))/2.d0
     &          -DC(n)*(TA(n)+TC(n))
      ENDIF
      IF(n.gt.2) THEN
         DO i=1,n
            PTZ(i)= DL(i-1)*PCZ(i)*(AL(i-1)-ALC(i-1))/2.d0 
     &          -DL(i)*(PCZ(i+1)-PCZ(i))*(AL(i)-ALC(i))/2.d0
     &          +DL(i+1)*PCZ(i+1)*(AL(i+1)-ALC(i+1))/2.d0
     &          -DC(i)*(TA(i)-TC(i))
         ENDDO
      ENDIF

      RETURN

      ENTRY ITANH3(N,Z0,Z,T0,TM,DL,C,CAV,DLC,ALC,DC
     & ,TC,PCZ)
C
C     CONVERT 'T' ARRAY INPUT TO 'C' ARRAY
c     Input: 
c            N+2 = number of profile points 
c            Z0 = the independent variable at the lowest profile point
c            Z = array of independent variables in profile
c            T0 = the dependent variable at the lowest profile point
c            TM = array of dependent variable in profile
c            DL = array of transition lengths in profile
c     Output:
c            C = array of slopes for profile
C            CAV = the average of the first and last values in the C array
C            DLC, ALC, DC, and TC are arrays that are based on the profiles 
C                 and are calculated once to avoid recalculating each time.
!            pCz = partial C(i)/partial z(i)

      TIM1=T0
      ZIM1=Z0
      NP1=N+1
C
      DO 10 I=1,NP1
         TI=TM(I)
         ZI=Z(I)
         C(I)=(TI-TIM1)/(ZI-ZIM1)
         pCz(I) = -C(I)/(ZI-ZIM1)
         TIM1=TI
10       ZIM1=ZI
C
      CAV = 0.5d0 * (C(1) + C(N + 1)) 
C
      DO 20 I = 1,N
      DLC(I) = DL(I) * (C(I + 1) - C(I)) / 2.d0 
      ALC(I) = ALCOSH((Z(I)-Z0) / DL(I))
      DC(I) = (C(I + 1) - C(I)) / 2.d0 
      TC(I) = dTANH((Z(I)-Z0) / DL(I))
20    CONTINUE

      RETURN

C ----------------------------------------------------------------------!-------
      END SUBROUTINE FTANH3
C ----------------------------------------------------------------------!-------




C ----------------------------------------------------------------------!-------
C from S4
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE GAMANG(THA,PHA)
C ----------------------------------------------------------------------!-------
C CALCULATE THE SIDE OF A SPHERICAL TRIANGLE BY LAW OF COSINES
C ----------------------------------------------------------------------!-------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/SPHGAM/SINLM0,COSLM0,GPH0,COSPHD,SINTH,COSTH
      COMMON/SPHGAM/GAMFUN,PGMTH,PGMPH,PGMTHTH,PGMPHPH,PGMTHPH
      COSPHD=dCOS(PHA-GPH0)
      SINTH=dSIN(THA)
      COSTH=dCOS(THA)
      GAMFUN=SINLM0*COSTH+COSLM0*SINTH*COSPHD
      PGMTH=-SINLM0*SINTH+COSLM0*COSTH*COSPHD
      PGMPH=-COSLM0*SINTH*dSIN(PHA-GPH0)
      PGMTHTH=-(SINLM0*COSTH+COSLM0*SINTH*COSPHD)
      PGMPHPH=-COSLM0*SINTH*COSPHD
      PGMTHPH=-COSLM0*COSTH*dSIN(PHA-GPH0)
      RETURN
C ----------------------------------------------------------------------!-------
       END SUBROUTINE GAMANG
C ----------------------------------------------------------------------!-------



C ----------------------------------------------------------------------!-------
C later additions
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE PLTNOTE(WLABEL,IDSU)                            ! 15 Mar 2010
C ----------------------------------------------------------------------!-------
C Subroutine to add labels for contour plots

      implicit double precision (a-h,o-z)

      COMMON /DD/ INT,IOR,IT,IS,IC,ICC,IX,IY

      INTEGER IDSU(20)

      DOUBLE PRECISION WLABEL(24)

      CHARACTER NUMSTG*80

      CHARACTER*80 LABEL,CHID,id
      common/ident/id

      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      WRITE(LABEL,'(a)') ID

      label = NUMSTG(idsu,20,'(20a4)')
c     print *, ' LABEL = ',LABEL

      IOR=0
      INT=1
      IT=0
      IS=3

      IS=1   
      IS=0   
      INT=0

      CHID=LABEL

      LABEL=CHID(1:10)
      IX=WLABEL( 1)
      IY=WLABEL( 2)
      IF(IY.NE.0) CALL DDTEXT (2,LABEL)

      LABEL=CHID(11:20)
      IX=WLABEL( 3)
      IY=WLABEL( 4)
      IF(IY.NE.0) CALL DDTEXT (2,LABEL) 

      LABEL=CHID(21:30)
      IX=WLABEL( 5)
      IY=WLABEL( 6)
      IF(IY.NE.0) CALL DDTEXT (2,LABEL) 

      LABEL=CHID(31:40)
      IX=WLABEL( 7)
      IY=WLABEL( 8)
      IF(IY.NE.0) CALL DDTEXT (2,LABEL) 

      LABEL=CHID(41:50)
      IX=WLABEL( 9)
      IY=WLABEL(10)
      IF(IY.NE.0) CALL DDTEXT (2,LABEL) 

      LABEL=CHID(51:60)
      IX=WLABEL(11)
      IY=WLABEL(12)
      IF(IY.NE.0) CALL DDTEXT (2,LABEL)

      LABEL=CHID(61:70)
      IX=WLABEL(13)
      IY=WLABEL(14)
      IF(IY.NE.0) CALL DDTEXT (2,LABEL) 

      LABEL=CHID(71:80)
      IX=WLABEL(15)
      IY=WLABEL(16)
      IF(IY.NE.0) CALL DDTEXT (2,LABEL) 

      LABEL=CHID(31:40)
      IX=WLABEL(17)
      IY=WLABEL(18)
C      IF(IY.NE.0) CALL DDTEXT (2,LABEL) 

      LABEL=CHID(41:50)
      IX=WLABEL(19)
      IY=WLABEL(20)
C      IF(IY.NE.0) CALL DDTEXT (2,LABEL) 

      LABEL=CHID(51:60)

      IF(IY.NE.0) CALL DDTEXT (2,LABEL) 

      IS=1
      INT=0
      RETURN
C ----------------------------------------------------------------------!-------
       END SUBROUTINE PLTNOTE
C ----------------------------------------------------------------------!-------



C     ***QUARTIC************************************************25.03.98
C     Solution of a quartic equation
C     ref.: J. E. Hacke, Amer. Math. Monthly, Vol. 48, 327-328, (1941)
C     NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
C     ------------------------------------------------------------------!-------
C     ******************************************************************
C     dd(0:4)     (i)  vector containing the polynomial coefficients
C     sol(1:4)    (o)  results, real part
C     soli(1:4)   (o)  results, imaginary part
C     Nsol        (o)  number of real solutions 
C     ==================================================================
      subroutine quartic(dd,sol,soli,Nsol)
C     ==================================================================
      implicit double precision (a-h,o-z)
      dimension dd(0:4),sol(4),soli(4)
      dimension AA(0:3),z(3)
C
      Nsol = 0
      a = dd(4)
      b = dd(3)
      c = dd(2)
      d = dd(1)
      e = dd(0)
C
      if (dd(4).eq.0.d+0) then
	write(6,*)'ERROR: NOT A QUARTIC EQUATION'
	return
      endif
C
      p = (-3.d+0*b**2 + 8.d+0*a*c)/(8.d+0*a**2)
      q = (b**3 - 4.d+0*a*b*c + 8.d+0*d*a**2)/(8.d+0*a**3)
      r = (-3.d+0*b**4 + 16.d+0*a*b**2*c - 64.d+0*a**2*b*d + 
     &      256.d+0*a**3*e)/(256.d+0*a**4)
C
C     solve cubic resolvent
      AA(3) =  8.d+0
      AA(2) = -4.d+0*p 
      AA(1) = -8.d+0*r
      AA(0) =  4.d+0*p*r - q**2
      call cubic(AA,z,ncube)
C      
      zsol = -1.d+99
      do 5 i=1,ncube
 5      zsol = max(zsol,z(i))
      z(1) = zsol
      xK2 = 2.d+0 * z(1) - p
      xK  = sqrt(xK2)
C ----------------------------------------------------------------------!-------
      if (xK.eq.0.d+0) then
        xL2 = z(1)**2 - r
	if (xL2.lt.0.d+0) then
	  write(6,*)'Sorry, no solution'
	  return
        endif
	xL  = sqrt(xL2)
      else
        xL = q/(2.d+0 * xK)
      endif
C ----------------------------------------------------------------------!-------
      sqp = xK2 - 4.d+0*(z(1) + xL)
      sqm = xK2 - 4.d+0*(z(1) - xL)
C
      do 10 i=1,4
 10     soli(i) = 0.d+0
      if       (sqp.ge.0.d+0 .and. sqm.ge.0.d+0) then
	sol(1) = 0.5d+0*( xK + sqrt(sqp))
	sol(2) = 0.5d+0*( xK - sqrt(sqp))
	sol(3) = 0.5d+0*(-xK + sqrt(sqm))
	sol(4) = 0.5d+0*(-xK - sqrt(sqm))
	Nsol = 4
      else if  (sqp.ge.0.d+0 .and. sqm.lt.0.d+0) then
	sol(1) =  0.5d+0*(xK + sqrt(sqp))
	sol(2) =  0.5d+0*(xK - sqrt(sqp))
	sol(3) = -0.5d+0*xK 
	sol(4) = -0.5d+0*xK 
	soli(3) =  sqrt(-0.25d+0 * sqm)
	soli(4) = -sqrt(-0.25d+0 * sqm)
	Nsol = 2
      else if  (sqp.lt.0.d+0 .and. sqm.ge.0.d+0) then
	sol(1) = 0.5d+0*(-xK + sqrt(sqm))
	sol(2) = 0.5d+0*(-xK - sqrt(sqm))
	sol(3) =  0.5d+0*xK 
	sol(4) =  0.5d+0*xK 
	soli(3) =  sqrt(-0.25d+0 * sqp)
	soli(4) = -sqrt(-0.25d+0 * sqp)
	Nsol = 2
      else if  (sqp.lt.0.d+0 .and. sqm.lt.0.d+0) then
	sol(1) = -0.5d+0*xK 
	sol(2) = -0.5d+0*xK 
	soli(1) =  sqrt(-0.25d+0 * sqm)
	soli(2) = -sqrt(-0.25d+0 * sqm)
	sol(3) =  0.5d+0*xK 
	sol(4) =  0.5d+0*xK 
	soli(3) =  sqrt(-0.25d+0 * sqp)
	soli(4) = -sqrt(-0.25d+0 * sqp)
	Nsol = 0
      endif
      do 20 i=1,4
 20     sol(i) = sol(i) - b/(4.d+0*a)
C
      return
C ----------------------------------------------------------------------!-------
      END SUBROUTINE QUARTIC
C ----------------------------------------------------------------------!-------

C     ***CUBIC************************************************08.11.1986
C     Solution of a cubic equation
C     Equations of lesser degree are solved by the appropriate formulas.
C     The solutions are arranged in ascending order.
C     NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING
C     ******************************************************************
C     A(0:3)      (i)  vector containing the polynomial coefficients
C     X(1:L)      (o)  results
C     L           (o)  number of valid solutions (beginning with X(1))
C ----------------------------------------------------------------------!-------
      SUBROUTINE CUBIC(A,X,L)
C ----------------------------------------------------------------------!-------
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(0:3),X(3),U(3)
      PARAMETER(PI=3.1415926535897932D+0,THIRD=1.D+0/3.D+0)
      INTRINSIC MIN,MAX,ACOS
C
C     define cubic root as statement function
      CBRT(Z)=SIGN(ABS(Z)**THIRD,Z)
C
C     ==== determine the degree of the polynomial ====
C
      IF (A(3).NE.0.D+0) THEN
C
C       cubic problem
	W=A(2)/A(3)*THIRD
	P=(A(1)/A(3)*THIRD-W**2)**3
	Q=-.5D+0*(2.D+0*W**3-(A(1)*W-A(0))/A(3))
	DIS=Q**2+P
	IF (DIS.LT.0.D+0) THEN
C         three real solutions!
C         Confine the argument of ACOS to the interval [-1;1]!
	  PHI=ACOS(MIN(1.D+0,MAX(-1.D+0,Q/SQRT(-P))))
	  P=2.D+0*(-P)**(5.D-1*THIRD)
	  DO 100 I=1,3
  100       U(I)=P*COS((PHI+DBLE(2*I)*PI)*THIRD)-W
	  X(1)=MIN(U(1),U(2),U(3))
	  X(2)=MAX(MIN(U(1),U(2)),MIN(U(1),U(3)),MIN(U(2),U(3)))
	  X(3)=MAX(U(1),U(2),U(3))
	  L=3
	ELSE
C         only one real solution!
	  DIS=SQRT(DIS)
	  X(1)=CBRT(Q+DIS)+CBRT(Q-DIS)-W
	  L=1
	END IF
C
      ELSE IF (A(2).NE.0.D+0) THEN
C
C       quadratic problem
	P=5.D-1*A(1)/A(2)
	DIS=P**2-A(0)/A(2)
	IF (DIS.GE.0.D+0) THEN
C         two real solutions!
	  X(1)=-P-SQRT(DIS)
	  X(2)=-P+SQRT(DIS)
	  L=2
	ELSE
C         no real solution!
	  L=0
	END IF
C
      ELSE IF (A(1).NE.0.D+0) THEN
C
C       linear equation
	X(1)=-A(0)/A(1)
	L=1
C
      ELSE
C       no equation
	L=0
      END IF
C
C     ==== perform one step of a newton iteration in order to minimize
C          round-off errors ====
      DO 110 I=1,L
	X(I)=X(I)-(A(0)+X(I)*(A(1)+X(I)*(A(2)+X(I)*A(3))))
     *  /(A(1)+X(I)*(2.D+0*A(2)+X(I)*3.D+0*A(3)))
  110 CONTINUE
      RETURN
C ----------------------------------------------------------------------!-------
      END SUBROUTINE CUBIC
C ----------------------------------------------------------------------!-------



