C file S7.FOR
C ----------------------------------------------------------------------!-------
C 30 Apr 2009 - replaced REAL*8 with DOUBLE PRECISION
C 10 Sep 2008 - Made data types consistent.
C 13 Mar 2008 - Corrected length in common block /PCONST/
C 13 Mar 2008 - Brought common block /HDR/ up to date
C 24 Mar 2006 - Moved COPYVL, CROSSP, CTRIM, SETUP, IZPZ, ZPZ, GCIRCLE, MATRIXM,
C               NORMV from contour1
C 24 Mar 2006 - Moved SWAP2, BINARY, FUD from eigen2
C 24 Mar 2006 - Moved XLOR from eigen
C 24 Mar 2006 - Moved C from Eigen3 and pcprof1
C ----------------------------------------------------------------------!-------
CONTAINS FUNCTION C, FUNCTION XLOR, FUNCTION BINARY, FUNCTION FUD, 
CONTAINS SUBROUTINE SWAP2, SUBROUTINE COPYVL, SUBROUTINE CROSSP, FUNCTION CTRIM, 
CONTAINS SUBROUTINE SETUP, SUBROUTINE IZPZ, SUBROUTINE ZPZ, FUNCTION GCIRCLE, 
CONTAINS SUBROUTINE MATRIXM, SUBROUTINE NORMV 
C ----------------------------------------------------------------------!-------
      FUNCTION C(T)
C ----------------------------------------------------------------------!-------
C     TRUNCATE INPUT STRING TO 81 CHARACTERS
C     ENTRY CD RETURNS STRING WITH '$' APPENDED AFTER LAST NON-BLANK
C     CHARACTER TO CONFORM TO DISSPLA PACKAGE REQUIREMENTS
C
      PARAMETER (NCHRS=81)
      CHARACTER C*(NCHRS),CD*(NCHRS),T*(*)
      INTEGER STRIM
      C=T
      RETURN
C
C ----------------------------------------------------------------------!-------
      ENTRY CD(T)
C ----------------------------------------------------------------------!-------
C     BE SURE STRING IS TERMINATED BY A $
      C=T
C     LIMIT THE NUMBER OF CHARACTERS TO PARAMETER 'NCHRS'
      NC=MIN0(NCHRS,STRIM(C))
      NCP=MIN0(NCHRS,NC+1)
      IF(C(NC:NC).NE.'$') C(NCP:NCP)='$'
C ----------------------------------------------------------------------!-------
      END FUNCTION C
C ----------------------------------------------------------------------!-------
C			Function XLOR
C ----------------------------------------------------------------------!-------
C     EXCLUSIVE OR 
 
      LOGICAL FUNCTION XLOR(XL1,XL2)
      LOGICAL XL1,XL2
C
      XLOR=(XL1.AND..NOT.XL2).OR.(.NOT.XL1.AND.XL2)
C ----------------------------------------------------------------------!-------
      END FUNCTION XLOR
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
C			Function BINARY
C VALUE OF FUNCTION IS 1 IF ARGUMENT GREATER OR EQUAL TO ZERO, -1 OTHERWISE
C ----------------------------------------------------------------------!-------
 
       double precision FUNCTION BINARY(X)
	
 	DOUBLE PRECISION x
 
       BINARY=1.0
       IF(X.GE.0.0) RETURN
       BINARY=-1.0
C ----------------------------------------------------------------------!-------
       END FUNCTION BINARY
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
C			Function FUD
C
C ----------------------------------------------------------------------!-------
C 13 September 2005 - RETURNS VALUE SLIGHTLY HIGH IN ALL CASES FOR
C                     CONVERTING TO INTEGERS TO AVOID TRUNCATION PROBLEM

      FUNCTION FUD(X)
 
      implicit double precision (a-h,o-z)
 
C     RETURN VALUE SLIGHTLY LOW FOR NEGATIVE VALUES, SLIGHTLY HIGH
C     FOR POSITIVE VALUES. USED WHEN ROUNDOFF PROBLEMS MESS COMPARES.
C
      FUD=0.0
C***      IF(X.LT.0.0) FUD=1.001*X
C***      IF(X.GT.0.0) FUD=.999*X     !   <-- NO!
C*** Changed by A Newhall    Bug?
      FUD=1.001*X
 
      RETURN
C ----------------------------------------------------------------------!-------
      END FUNCTION FUD
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
C			Subroutine SWAP2
C ----------------------------------------------------------------------!-------
C SWAPS TWO PAIRS OF NUMBERS 
C ----------------------------------------------------------------------!-------
      SUBROUTINE SWAP2(X1,X2,Y1,Y2)
C ----------------------------------------------------------------------!-------
 
	implicit double precision (a-h,o-z)
 
C
        T=X1
        X1=X2
        X2=T
        T=Y1
        Y1=Y2
        Y2=T
C ----------------------------------------------------------------------!-------
      END SUBROUTINE SWAP2
C ----------------------------------------------------------------------!-------
 
C ----------------------------------------------------------------------!-------
      SUBROUTINE COPYVL(V,W)
C ----------------------------------------------------------------------!-------
C COPIES A VECTOR ONTO A COLUMN OF A MATRIX.
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

      DOUBLE PRECISION V(3),W(3,3)
C
      DO 10 I=1,3
10    W(1,I)=V(I)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE COPYVL
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE CROSSP(V,W,R)
C ----------------------------------------------------------------------!-------
C
C     COMPUTE CROSS PRODUCT OF VECTORS V AND W , RESULT IN R
C
      implicit double precision (a-h,o-z)

      DOUBLE PRECISION V(3),W(3),R(3)
C
      R(1)=V(2)*W(3)-V(3)*W(2)
      R(2)=-(V(1)*W(3)-V(3)*W(1))
      R(3)=V(1)*W(2)-V(2)*W(1)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE CROSSP
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      INTEGER FUNCTION CTRIM(C)
C ----------------------------------------------------------------------!-------
C FINDS LENGTH OF CHARACTER STRING EXCLUDING TRAILING BLANKS. 
C ----------------------------------------------------------------------!-------
      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
      CTRIM=1
      RETURN
200   CTRIM=I
C ----------------------------------------------------------------------!-------
      END FUNCTION CTRIM
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE SETUP
C ----------------------------------------------------------------------!-------
C     OBTAIN DATE-TIME, READ USER IDENTIFICATION AND COMPUTE CONSTANTS.
      implicit double precision (a-h,o-z)

C     COMMON DECK "RAYDEV" INSERTED HERE
C     DEVICE ASSIGNED TO RAYTRC INPUT FILE
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVTXT,NDEVBIN
C
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 "HDR" INSERTED HERE
      CHARACTER*10 INITID*80,DAT,TOD
      COMMON/HDR/SEC,stsec                            !Rev 7/13/90 ! 13 Mar 2008
      COMMON/HDRC/INITID,DAT,TOD
      PARAMETER (NWARSZ=1000)
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
C
ccccc      CALL SYSDAT(DAT)
ccccc      CALL SYSTIM(TOD)
ccccc      CALL SYSSEC(SEC)
C     READ IN USERS NAME AND TELEPHONE EXTENSION FOR IDENTIFYING MICROF
C     PLOTS.
ccccc      READ(NRYIND,1100,END=10000) INITID
ccccc1100  FORMAT(A)
ccccc10000 CALL DDINIT(8,INITID)

CONPLT.21
C********* INITIALIZE MATH CONSTANTS AND THE W ARRAY
C      CALL SETW                            ! 2 Mar 2006
      gronk = setw(xx)                      ! 2 Mar 2006
C ----------------------------------------------------------------------!-------
      END SUBROUTINE SETUP
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE IZPZ(THETL,PHIL,THETR,PHIR)
C ----------------------------------------------------------------------!-------
C INITIALIZATION FOR ZPZ 
C ----------------------------------------------------------------------!-------

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

      DOUBLE PRECISION KPX,KPY,KPZ,IPX,IPY,IPZ
      DOUBLE PRECISION LX,LY,LZ
      DOUBLE PRECISION L(3),R(3),KP(3),IP(3),V(3),RS(3)
C
      COMMON/ZPZC/LX,LY,LZ,RX,RY,RZ,KPX,KPY,KPZ,IPX,IPY,IPZ
      COMMON/ZPZC/VX,VY,VZ,RSX,RSY,RSZ
      COMMON/ZPZC/BETA,A(3,3)
C
      EQUIVALENCE (L,LX),(R,RX),(KP,KPX),(IP,IPX),(V,VX)
      EQUIVALENCE (RS,RSX)
C
C
      S=SIN(THETL)
      LX=S*COS(PHIL)
      LY=S*SIN(PHIL)
      LZ=COS(THETL)
C
      S=SIN(THETR)
      RX=S*COS(PHIR)
      RY=S*SIN(PHIR)
      RZ=COS(THETR)
C
C     DETERMINE PRIMED Z AXIS BY CROSS PRODUCT OF L AND R
      CALL CROSSP(L,R,KP)
      CALL NORMV(KP)
      IF(ABS(ABS(KPZ)-1.0).GT.1.E-10) GO TO 10
C
C     ONLY A SIMPLE Z AXIS ROTATION IS REQUIRED, SKIP MATRIX COMPUTATION
C     MAKE IT A CLEAN 90 DEGREES, MINUS ANGLES FOR SOUTH POLE CASE
C
      KPZ=SIGN(1.0d0,KPZ)                                          ! 10 Sep 2008
      BETA=PHIL*KPZ
      RETURN
C
10    CALL CROSSP(KP,L,IP)
      CALL COPYVL(L,A)
      CALL COPYVL(IP,A(2,1))
      CALL COPYVL(KP,A(3,1))
C
C ----------------------------------------------------------------------!-------
      END SUBROUTINE IZPZ
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE ZPZ(THETAP,PHIP,THETA,PHI)
C ----------------------------------------------------------------------!-------
C
C     'ZPZ' CONVERTS THE SPHERICAL COORDINATES THETAP,PHIP
C     OF A POINT IN PRIMED SPACE TO COORDINATES THETA,PHI IN
C     UNPRIMED SPACE WHERE PRIMED SPACE IS LOCATED RELATIVE TO
C     UNPRIMED SPACE SUCH THAT POINTS THETL,PHIL AND THETR,PHIR
C     LIE ON THE EQUATOR. SEE IZPZ FOR INITIALIZATION. METHOD
C     USED IS TO CONVERT TO CARTESIAN COORDINATES, ROTATE TO UNPRIMED
C     SYSTEM AND CONVERT BACK TO SPHERICAL COORDINATES.
C
      implicit double precision (a-h,o-z)
C
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
      DOUBLE PRECISION KPX,KPY,KPZ,IPX,IPY,IPZ
      DOUBLE PRECISION LX,LY,LZ
C
      COMMON/ZPZC/LX,LY,LZ,RX,RY,RZ,KPX,KPY,KPZ,IPX,IPY,IPZ
      COMMON/ZPZC/VX,VY,VZ,RSX,RSY,RSZ
      COMMON/ZPZC/BETA,A(3,3)
C
      IF(ABS(KPZ).NE.1.0) GO TO 100
      THETA=THETAP
      IF(KPZ.LT.0.0) THETA=PI-THETAP
      PHI=PHIP+BETA
      RETURN
C
100   S=1.0
      IF(THETAP.NE.PID2) S=SIN(THETAP)
      VX=S*COS(PHIP)
      VY=S*SIN(PHIP)
      VZ=0.0
      IF(THETAP.NE.PID2) VZ=COS(THETAP)
C
      CALL MATRIXM(A,VX,RSX)
      THETA=ACOS(RSZ)
      PHI=ATAN2(RSY,RSX)
      IF(ABS(RSZ).EQ.1.0) PHI=0.0
C ----------------------------------------------------------------------!-------
      END SUBROUTINE ZPZ
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      FUNCTION GCIRCLE(LLAT,LLON,RLAT,RLON)
C ----------------------------------------------------------------------!-------
C GREAT CIRCLE DISTANCE 
C ----------------------------------------------------------------------!-------
      implicit double precision (a-h,o-z)

      DOUBLE PRECISION LLAT,LLON
C
      CTHR=SIN(RLAT)
      STHR=COS(RLAT)
      CTHL=SIN(LLAT)
      STHL=COS(LLAT)
      CLR=CTHL*CTHR+STHL*STHR*COS(RLON-LLON)
      GCIRCLE=ACOS(CLR)
C ----------------------------------------------------------------------!-------
      END FUNCTION GCIRCLE
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE MATRIXM(A,V,R)
C ----------------------------------------------------------------------!-------
C MATRIX MULTPLICATION 
C ----------------------------------------------------------------------!-------
      DOUBLE PRECISION A(3,3),V(3),R(3)
C
      DO 10 I=1,3
      R(I)=0.0
      DO 10 J=1,3
10    R(I)=R(I)+A(J,I)*V(J)
C ----------------------------------------------------------------------!-------
      END SUBROUTINE MATRIXM
C ----------------------------------------------------------------------!-------
C ----------------------------------------------------------------------!-------
      SUBROUTINE NORMV(V)
C ----------------------------------------------------------------------!-------
C
C     COMPUTE NORMLIZED VECTOR COMPONENTS FOR VECTOR V
C
      DOUBLE PRECISION V(3)
C
      R=V(1)*V(1)+V(2)*V(2)+V(3)*V(3)
      IF(R.LT.1.E-10) RETURN
      R=1.0/SQRT(R)
      DO 10 I=1,3
10    V(I)=R*V(I)
C
C ----------------------------------------------------------------------!-------
      END SUBROUTINE NORMV
C ----------------------------------------------------------------------!-------
