      SUBROUTINE TROPLV2 (PR, HEIGHT, TEMPS, N, TROTMP, TROPRS)
      IMPLICIT NONE
C
C***********************************************************************
C
C   NAME- TROPLV
C
C   LANGUAGE- VS FORTRAN      TYPE- SUBROUTINE      COMPUTER- NAS9000
C
C   VERSION- 1.0      DATE- 2/27/92      PROGRAMMER-  DAVE MOREL (SMSRC)
C   VERSION- 2.0      DATE- 3/22/95      PROGRAMMER-  CARLA NEIGH(SMSRC)
C    ** ATOVS-68 : CHECKS ADDED FOR MISSING VALUES OF HEIGHT AND LAPSE
C                  RATE.  ASSUMPTION ADDED THAT MISSING VALUES CAN
C                  OCCUR, BUT ONLY AT THE TOP AND BOTTOM OF A PROFILE.
C                  ASSUMPTION FOR 40 TOVS LEVELS CHANGED TO 42 LEVELS.
C   VERSION- 3.0    DATE- 05/14/96   PROGRAMMER- ALEX POZDNYAK, FDC
C      FOR REHOSTING TO THE CRAY, CHANGED ALL I*4 DECLARATIONS TO
C      INTEGER, ALL L*4 TO LOGICAL, AND ALL R*4 DECLARATIONS TO REAL.
C   VERSION- 4.0    DATE- 02/27/03   PROGRAMMER- Rickey Rollins (CSC)
C        Replace R4SPV equivalence test with ISR4SPV function 
C
C   FUNCTION- THIS SUBROUTINE DETERMINES THE TROPOPAUSE LEVEL AND
C   RETURNS ITS TEMPERATURE (K) AND PRESSURE (MB).
C
C   CALLING SEQUENCE PARAMETERS-
C
C   NAME    TYPE  I/O    DESCRIPTION
C   ----    ----  ---    -----------
C   PR      REAL   I     VECTOR OF N LEVEL PRESSURES (MB),ASSUMED TO BE
C                          THE 42 TOVS LEVELS.
C   HEIGHT  REAL   I     VECTOR OF N GEOMETRIC HEIGHT IN (KM),
C                          CORRESPONDING TO THE LEVELS IN VECTOR PR.
C   TEMPS   REAL   I     VECTOR OF N TEMPERATURES (K),
C                          CORRESPONDING TO THE LEVELS IN VECTOR PR.
C   N       INT    I     NUMBER OF PRESSURE LEVELS, USED ONLY TO DEFINE
C                          THE DIMENSION OF THE PR, HEIGHT AND TEMPS
C                          VECTORS.
C   TROTMP  REAL   O     TROPOPAUSE TEMPERATURE (K).
C   TROPRS  REAL   O     TROPOPAUSE PRESSURE (MB).
C
      INTEGER     N
      REAL        PR(N), HEIGHT(N), TEMPS(N), TROTMP, TROPRS
C
C   COMMON AREAS-  NONE
C
C   LOCAL VARIABLES
C   NAME      TYPE  DESCRIPTION
C   ----      ----  -----------
C   I,J,K      INT  DO LOOP INDEX
C   L          INT  HOLD VARIABLE FOR INDEX J-1
C   LPSRTE    REAL  VECTOR OF MAX 100 LAPSE RATES
C   LPSRAV    REAL  STATEMENT FUNCTION TO COMPUTE AVERAGE LAPSE RATE
C   OKFLAG     LOG  FLAG SET TO TRUE IF LAPSE RATE <= 2 K/KM
C   AVLAPR    REAL  AVERAGE LAPSE RATE OVER 2 KM
C   DTAU      REAL  THICKNESS IN KM OF THE HEIGHT SUBINTERVAL CONTAINING
C                   THE TROPOPAUSE
C   DIVDF1    REAL  DIVIDED DIFFERENCE USED TO COMPUTE POLYNOMIAL COEFFS
C   DIVDF1    REAL  DIVIDED DIFFERENCE USED TO COMPUTE POLYNOMIAL COEFFS
C   C1        REAL  CUBIC HERMITE POLYNOMIAL COEFFICIENT
C   C2        REAL  CUBIC HERMITE POLYNOMIAL COEFFICIENT
C   C3        REAL  CUBIC HERMITE POLYNOMIAL COEFFICIENT
C   C4        REAL  CUBIC HERMITE POLYNOMIAL COEFFICIENT
C   A         REAL  QUADRATIC EQUATION COEFFICIENT
C   B         REAL  QUADRATIC EQUATION COEFFICIENT
C   C         REAL  QUADRATIC EQUATION COEFFICIENT
C   DISCR     REAL  DISCRIMINANT OF SOLUTION TO QUADRATIC EQUATION
C   Q         REAL  INTERMEDIATE QUANTITY IN SOLVING QUADRATIC EQUATION
C   X1        REAL  1ST ROOT OF QUADRATIC EQUATION
C   X2        REAL  2ND ROOT OF QUADRATIC EQUATION
C   DELZ      REAL  VERTICAL DISTANCE BETWEEN THE REFINED VALUE OF THE
C                   TROPOPAUSE AND THE HEIGHT OF THE UPPER BOUND OF THE
C                   SUBINTERVAL.
C   CONST     REAL  COMPILER CONSTANT, VALUE OF GROUPED CONSTANTS IN THE
C                   SOLUTION TO THE HYPSOMETRIC EQUATION:
C                         (2*G/R) * 1000(M/KM)
C
C   SUBPROGRAMS CALLED-  LAPSER, LAPRTE
C
C   FUNCTIONS USED-  EXP, MIN, SIGN, SQRT
C
C   EXIT STATES- NONE
C
C   ASSUMPTIONS - 1) THE INPUT DATA ARE AT THE 42 TOVS LEVELS, IN ORDER
C                    OF THE HIGHEST LEVEL 1ST AND 1000 MB AT LEVEL 40
C                    (LEVELS 41 AND 42 ARE NOT USED).
C                 2) THE TEMPERATURES AND HEIGHTS CAN BE MISSING VALUES,
C                    HOWEVER IT IS ASSUMED THAT THESE VALUES OCCUR ONLY
C                    AT THE TOP AND/OR BOTTOM OF THE PROFILE.  THAT IS,
C                    BETWEEN THE HIGHEST AND LOWEST NON-MISSING VALUES,
C                    THERE ARE NO MISSING VALUES.  IT IS ASSUMED THAT
C                    FOR ANY LEVEL, HEIGHT AND TEMPERATURE ARE EITHER
C                    BOTH MISSING OR BOTH PRESENT.
C***********************************************************************
C
      INTEGER    I, J, K, L
      REAL       LPSRTE(100),AVLAPR,CONST, LPSRAV
      PARAMETER (CONST = 2.*1000.*9.8/287.04)
      REAL       DTAU,DIVDF1,DIVDF3,C1,C2,C3,C4,A,B,C,DISCR,Q,X1,X2,DELZ
      LOGICAL    OKFLAG
      LOGICAL   ISR4SPV
C
      INCLUDE 'SPVAL'
C
C  STATEMENT FUNCTION TO COMPUTE AVERAGE LAPSE RATE
C
      LPSRAV(J,I) = -(TEMPS(J) - TEMPS(I))/(HEIGHT(J) - HEIGHT(I))
C
      TROTMP  = R4SPV
      TROPRS  = R4SPV
C
C  THE HEIGHTS ARE ASSUMED TO BE IN GEOMETRIC KILOMETERS, THE LAPSE
C  RATE IS DEGREES PER KM
C
C  CALCULATE THE LAPSE RATES AT EACH LEVEL
C
      CALL LAPSER(HEIGHT,TEMPS,N,LPSRTE)
ccft  CALL LAPRTE(HEIGHT,TEMPS,N,LPSRTE)
C
C  TO OBTAIN THE PROPER TROPOPAUSE LEVEL TO CHOOSE THE TEMPERATURE AND
C  PRESSURE, WE MUST DETERMINE THE LEVEL BY USING THE WMO DEFINITION
C  (E.G., PAGE B5-27 OF THE FEDERAL METEOROLOGICAL HANDBOOK #3, DATED
C   1-1-69.)
C
C  OBTAIN ›I:I+1| SUBINTERVAL OF PROFILE  CONTAINING THE TROPOPAUSE
C
C     LOOP FROM 500 MB TO 30 MB, WHICH CORRESPOND TO TOVS PRESSURE LEVEL
C     31 AND 15, RESPECTIVELY
C
      DO 30 I = 31, 15, -1
C
         IF (ISR4SPV(LPSRTE(I))) GOTO 40
C
C        FIND THE LOWEST LEVEL AT WHICH THE LAPSE RATE DECREASES TO
C        2 K/KM
C
         IF(LPSRTE(I) .LE. 2.0) THEN
c        IF(LPSRTE(I) .LE. 3.0) THEN
C
C           THE AVERAGE LAPSE RATE FROM THIS LEVEL TO ANY POINT
C           WITHIN THE NEXT HIGHER TWO KM MUST NOT EXCEED 2 K/KM
C
            DO 10 J = I-2, 14, -1
               IF (HEIGHT(J) .EQ. R4SPV) GOTO 30
               IF((HEIGHT(J)-HEIGHT(I))  .GT. 2.0) GO TO 20
   10       CONTINUE
   20       AVLAPR = LPSRAV(J,I)
            IF(AVLAPR .LE. 2.0) THEN
C
C  CHECK IF CANDIDATE TROPOPAUSE UPPER BOUND IS 500 MB.  IF SO, TROPO-
C  PAUSE IS REALLY BELOW 500 MB
C
              IF(I.EQ.31) GO TO 40
C
C  TROPOPAUSE UPPER BOUND IS NOT 500 MB, TROPOPAUSE CRITERIA MET,
C  PROCEED TO LOCATE ACTUAL VALUE
C
              GO TO 220
             ENDIF
         ENDIF
   30 CONTINUE
C
C     THE TROPOPAUSE LEVEL WAS NOT FOUND BETWEEN 500MB & 30MB NEXT
C     LOOP FROM 500 MB TO 1000MB
C
   40 OKFLAG = .FALSE.
      DO 100 I = 31, 40
C
         IF (ISR4SPV(LPSRTE(I))) GOTO 200
C
C        FIND THE HIGHEST LEVEL AT WHICH THE LAPSE RATE DECREASES TO
C        2 K/KM
C
C        IF THE LAPSE RATE IS GREATER THAN 2 K/KM AND WAS PREVIOUSLY
C        LESS THAN OR EQUAL TO 2 K/KM, THEN THE PREVIOUS LEVEL IS A
C        CANDIDATE FOR BEING A TROPOPAUSE UPPER BOUND
C
         IF((LPSRTE(I) .GT. 2.0).AND.OKFLAG) THEN
C
C           THE AVERAGE LAPSE RATE WITHIN THE NEXT HIGHER 2 KM MUST NOT
C           EXCEED 2 K/KM
C
            DO 50 J = I-3, 26, -1
               IF (ISR4SPV(HEIGHT(J))) GOTO 100
               IF((HEIGHT(J) - HEIGHT(I-1))  .LE. 2.0) GO TO 60
   50       CONTINUE
   60       AVLAPR = LPSRAV(J,I-1)
            IF(AVLAPR .LE. 2.0) THEN
C
C              THE AVERAGE LAPSE RATE MUST NOT EXCEED 3 K/KM OVER MORE
C              THAN 1 KM UP TO AND INCLUDING 100 MB (LEVEL 20)
C
               DO 90 K = I-1,21,-1
                 DO 70 J = K-1, 20, -1
                   IF (ISR4SPV(HEIGHT(J)).OR.ISR4SPV(HEIGHT(K)))GOTO 100
                   IF((HEIGHT(J) - HEIGHT(K)) .GT. 1.) THEN
                      IF((J+1).LT.K) THEN
                         L = J - 1
                         GO TO 80
                      ELSE
                         L = J
                      ENDIF
                   ENDIF
   70            CONTINUE
   80            AVLAPR = LPSRAV(L,K)
                 IF(AVLAPR .GT. 3.0) THEN
                   OKFLAG = .FALSE.
                   GO TO 100
                 ENDIF
   90          CONTINUE
C
C  TROPOPAUSE CRITERIA MET, SET I TO UPPER BOUND ONCE OUTSIDE THE LOOP
C
              GO TO 210
            ELSE
             OKFLAG = .FALSE.
            ENDIF
C
C  SET OKFLAG IF THE LAPSE RATE IS ENCOUNTERED THAT IS <= 2 K/KM
C
         ELSEIF(LPSRTE(I) .LE. 2.0) THEN
           OKFLAG = .TRUE.
         ENDIF
  100 CONTINUE
C
C  NO TROPOPAUSE FOUND, RETURN MISSING VALUES
C
  200 RETURN
C
C  UPPER BOUND FOR TROPOPAUSE FOUND, SEARCH THE LAYER BELOW FOR POINT
C  WHERE LAPSE RATE FALLS TO 2 K/KM
C
  210 I = I - 1
  220 TROTMP  = TEMPS(I)
      TROPRS  = PR(I)
C
C  FIT CUBIC HERMITE POLYNOMIAL TO SUBINTERVAL ›I:I+1| AND
C  SOLVE FOR THE POINT(S) IN THE SUBINTERVAL WHERE THE LAPSE RATE IS
C  2 K PER KM.  THE POLYNOMIAL IS
C
C  TEMPS(Z) = C1 + C2*(Z - Z(I+1)) + C3*(Z - Z(I+1))**2
C                                          + C4*(Z - Z(I+1))**3
C
C   (RECALLING THAT THE HEIGHTS, TEMPERATURES, PRESSURES AND LAPSE RATES
C   ARE STORED IN THEIR RESPECTIVE VECTORS WITH THE .1 MB LEVEL 1ST AND
C   THE 1000 MB LEVEL LAST.)
C
C  THE COEFFICIENTS OF CUBIC HERMITE POLYNOMIAL ON THE SUBINTERVAL ›I:I+
C  ARE COMPLETELY DEFINED FROM THE VALUES AND THE DERIVATIVES AT THE
C  SUBINTERVAL BOUNDARIES.  SEE DE BOOR, C.,1978: A PRACTICAL GUIDE TO
C  SPLINES, SPRINGER-VERLAG, 392 PP.  EQUATIONS (3) ON P 50.
C
      DTAU = HEIGHT(I) - HEIGHT(I+1)
      DIVDF1 = (TEMPS(I) - TEMPS(I+1))/DTAU
      DIVDF3 = - LPSRTE(I) - LPSRTE(I+1) - 2.*DIVDF1
      C1 = TEMPS(I+1)
      C2 = - LPSRTE(I+1)
      C3 = (DIVDF1 + LPSRTE(I+1) - DIVDF3)/DTAU
      C4 = DIVDF3/DTAU/DTAU
C
C     THE LAPSE RATE IS THE DERIVATIVE OF THE CUBIC POLYNOMIAL, AND
C     THEREFORE IS A QUADRATIC EQUATION
C
C    LPSRTE(Z) = - D(TEMPS)/DZ
C              = - C2 - 2.*C3*(Z - Z(I+1)) - 3.*C4*(Z - Z(I+1))**2
C
C     WE SOLVE QUADRATIC EQUATION (A*X**2 + B*X + C = 0) FOR  A LAPSE RA
C     0F 2, WHERE X = Z - Z(I+1).
C
      A = 3.*C4
      B = 2.*C3
      C = C2 + 2.
      DISCR = B*B - 4.*A*C
C
C  IF DISCRIMINANT IS NEGATIVE, THEN NO REAL ROOTS, SO RETURN
C  ALSO RETURN IF THE ROOTS DO NOT LIE IN THE HEIGHT SUBINTERVAL
C  I:I+1|
C
      IF(DISCR.LT.0.) RETURN
      Q = -(B + SIGN(SQRT(DISCR),B))/2.
      IF(Q.EQ.0.0) RETURN
      X2 = C/Q
      IF(A.NE.0) THEN
        X1 = Q/A
        IF((X2.LE.0.).AND.(X1.LE.0.)) RETURN
        IF((X2.LE.0.).AND.(X1.GT.0.)) THEN
          DELZ = X1
         ELSEIF((X2.GT.0.).AND.(X1.LE.0.)) THEN
          DELZ = X2
         ELSE
          DELZ = MIN(X1,X2)
         ENDIF
       ELSE
         IF(X2.LE.0.) RETURN
         DELZ = X2
       ENDIF
      IF(DELZ.GT.DTAU) RETURN
C
C  CHOOSE THE HIGHEST POINT AS THE TROPOPAUSE, AND EVALUATE
C  TEMPERATURE FROM THE POLYNOMIAL, AND CALCULATE PRESSURE FROM THE
C  HYPSOMETRIC EQUATION (SEE ANY METEOROLOGICAL TEXT).
C
C
      TROTMP  = C1 + DELZ*(C2 + DELZ*(C3 + C4*DELZ))

c     IF (TROTMP.LT.TEMPS(I).AND.TROTMP.LT.TEMPS(I+1)) THEN
c        TROTMP = (TEMPS(I) + TEMPS(I+1)) / 2.0
c     END IF
ccft
c     TROTMP = TEMPS(I)
      TROPRS  = PR(I+1)/EXP(CONST*DELZ/(TROTMP+TEMPS(I+1)))
C
      RETURN
      END
