      SUBROUTINE LAPSER(HEIGHT,TEMP,N,LPSRTE)
C
C***********************************************************************
C
C   NAME- LAPSER
C
C   LANGUAGE- FORTRAN 77      TYPE- SUBROUTINE      COMPUTER- CRAY   
C
C   VERSION- 1.0      DATE- 2/27/92      PROGRAMMER-  DAVE MOREL (SMSRC)
C   VERSION- 2.0      DATE- 1/13/95      PROGRAMMER-  CARLA NEIGH(SMSRC)
C    ATOVS-50:  ADD CODE TO CHECK FOR MISSING VALUES OF HEIGHT AND
C               TEMPERATURE.
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 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 PRODUCES THE LAPSE RATES AT THE NUMBER OF
C   GEOMETRIC HEIGHTS SPECIFIED IN THE CALLING ARGUMENT N.
C
C   CALLING SEQUENCE PARAMETERS-
C
C   NAME    TYPE  I/O    DESCRIPTION
C   ----    ----  ---    -----------
C   HEIGHT  REAL   I     VECTOR OF N GEOMETRIC HEIGHTS (KM).
C   TEMP    REAL   I     VECTOR OF N TEMPERATURE VALUES (K).
C   N       INT    I     NUMBER OF PRESSURE LEVELS.
C   LPSRTE  REAL   O     VECTOR OF N TO HOLD THE LAPSE RATES.  LAPSE
C                        RATE FOR A LEVEL WILL BE A MISSING VALUE IF
C                        IT COULD NOT BE COMPUTED DUE TO MISSING HEIGHT
C                        OR TEMPERATURE.
C
      INTEGER     N
      REAL        HEIGHT(N), LPSRTE(N), TEMP(N)
C
C   COMMON AREA PARAMETERS-
      INCLUDE 'SPVAL'
C
C   LOCAL VARIABLES
C   NAME      TYPE  DESCRIPTION
C   ----      ----  -----------
C   TOP        INT  HIGHEST LEVEL CONTAINING NON-MISSING VALUES.
C   BOTTOM     INT  LOWEST LEVEL CONTAINING NON-MISSING VALUES.
C   I          INT  DO LOOP INDEX.
C   NUMER1    REAL  FIRST NUMERATOR IN THE EQUATION.
C   NUMER2    REAL  SECOND NUMERATOR IN THE EQUATION.
C   DENOM1    REAL  FIRST DENOMINATOR IN THE EQUATION.
C   DENOM2    REAL  SECOND DENOMINATOR IN THE EQUATION.
C   EQUAT     REAL  VECTOR OF THE THREE PARTS OF THE FORMULA.
C
C   SUBPROGRAMS CALLED-  NONE
C
C   EXIT STATES- NONE
C
C   ASSUMPTIONS- MISSING VALUES CAN EXIST IN THE HEIGHT OR TEMPERATURE
C                ARRAYS ONLY AT TOP OR BOTTOM LEVELS.  ALL VALUES
C                BETWEEN THE LOWEST AND HIGHEST NON-MISSING VALUES WILL
C                BE NON-MISSING.
C
C   COMMENTS- THE REFERENCE FOR THE DIFFERENTIATION FORMULA USED IS
C             ELEMENTARY THEORY AND APPLICATION OF NMBEICAL ANALYSIS
C             BY MOURSUND & DURIS COPYRIGHT 1967.
C
C***********************************************************************
C
      INTEGER  I, TOP, BOTTOM
      REAL     NUMER1, NUMER2, DENOM1, DENOM2, EQUAT(3)
      LOGICAL   ISR4SPV
C
C     INITIALIZE LAPSE RATE TO MISSING VALUES.
C
      DO 10 I = 1, N
         LPSRTE(I) = R4SPV
   10 CONTINUE
C
C     FIND THE HIGHEST AND LOWEST LEVELS OF NON-MISSING VALUES OF
C     TEMPERATURE AND HEIGHT.
C
      TOP = I4SPV
      BOTTOM = I4SPV
C
      DO 20 I = 1, N
         IF ((.NOT. ISR4SPV(HEIGHT(I))) .AND. (.NOT. ISR4SPV(TEMP(I)))) THEN
            TOP = I
            GOTO 30
         ENDIF
   20 CONTINUE
C
   30 IF (TOP .EQ. I4SPV) RETURN
C
      DO 40 I = N, 1, -1
         IF ((.NOT. ISR4SPV(HEIGHT(I))) .AND. (.NOT. ISR4SPV(TEMP(I)))) THEN
            BOTTOM = I
            GOTO 50
         ENDIF
   40 CONTINUE
C
C     THREE LEVELS OF NON-MISSING VALUES ARE NEEDED TO COMPUTE THE
C     LAPSE RATES.
C
   50 IF ((BOTTOM - TOP) .LT. 2) RETURN
C
C     THE FIRST LEVEL USES THE SECOND AND THIRD LEVELS TO DETERMINE
C     THE LAPSE RATE USING THE FOLLOWING FORMULA:
C
C     LPSRTE(TOP) =-((((-(2*(HEIGHT(TOP+1)-HEIGHT(TOP)) +
C                  (HEIGHT(TOP+2)-HEIGHT(TOP+1)))
C                  / ((HEIGHT(TOP+1)-HEIGHT(TOP)) *
C                  * (HEIGHT(TOP+2)-HEIGHT(TOP))))) * TEMP(TOP)
C
C                 + ((((HEIGHT(TOP+2)-HEIGHT(TOP))  /
C                  ((HEIGHT(TOP+1)-HEIGHT(TOP)) *
C                  (HEIGHT(TOP+2)-HEIGHT(TOP+1))))) * TEMP(TOP+1)
C
C                 - ((((HEIGHT(TOP+1)-HEIGHT(TOP))  /
C                  ((HEIGHT(TOP+2)-HEIGHT(TOP)) *
C                  (HEIGHT(TOP+2)-HEIGHT(TOP+1)))))) * TEMP(TOP+2)
C
C     THE ABOVE FORMULA IS IMPLEMENTED AS FOLLOWS:
C
      NUMER1 = 2*(HEIGHT(TOP+1) - HEIGHT(TOP))
      NUMER2 = HEIGHT(TOP+2) - HEIGHT(TOP+1)
      DENOM1 = HEIGHT(TOP+1) - HEIGHT(TOP)
      DENOM2 = HEIGHT(TOP+2) - HEIGHT(TOP)
      EQUAT(1) = (-(NUMER1 + NUMER2) / (DENOM1 * DENOM2)) * TEMP(TOP)
      NUMER1 = HEIGHT(TOP+2) - HEIGHT(TOP)
      DENOM1 = HEIGHT(TOP+1) - HEIGHT(TOP)
      DENOM2 = HEIGHT(TOP+2) - HEIGHT(TOP+1)
      EQUAT(2) = (NUMER1 / (DENOM1 * DENOM2)) * TEMP(TOP+1)
      NUMER1 = HEIGHT(TOP+1) - HEIGHT(TOP)
      DENOM1 = HEIGHT(TOP+2) - HEIGHT(TOP)
      DENOM2 = HEIGHT(TOP+2) - HEIGHT(TOP+1)
      EQUAT(3) = (NUMER1 / (DENOM1 * DENOM2)) * TEMP(TOP+2)
      LPSRTE(TOP) = -(EQUAT(1) + EQUAT(2) - EQUAT(3))
C
C     THE SECOND THROUGH THE NEXT TO THE LAST LEVEL USES THE BEFORE
C     AND AFTER LEVELS TO DETERMINE THE LAPSE RATE USING THE FOLLOWING:
C
C         LPSRTE(I)= -((((-(HEIGHT(I+1)-HEIGHT(I)) /
C                 ((HEIGHT(I)-HEIGHT(I-1))*(HEIGHT(I+1)-HEIGHT(I-1)))))
C                 * TEMP(I-1))
C
C                - (((((HEIGHT(I)-HEIGHT(I-1))-(HEIGHT(I+1)-HEIGHT(I)))
C                 /((HEIGHT(I)-HEIGHT(I-1))*(HEIGHT(I+1)-HEIGHT(I))))) *
C                 TEMP(I))
C
C                 + ((((HEIGHT(I)-HEIGHT(I-1))  /
C                 ((HEIGHT(I+1)-HEIGHT(I-1))*(HEIGHT(I+1)-HEIGHT(I))))))
C                 * TEMP(I+1))
C
C     WHICH IS IMPLEMENTED AS FOLLOWS:
C
      DO 60 I = TOP+1, BOTTOM-1
        NUMER1 = -(HEIGHT(I+1) - HEIGHT(I))
        DENOM1 = HEIGHT(I) - HEIGHT(I-1)
        DENOM2 = HEIGHT(I+1) - HEIGHT(I-1)
        EQUAT(1) = (NUMER1 / (DENOM1 * DENOM2)) * TEMP(I-1)
        NUMER1 = HEIGHT(I) - HEIGHT(I-1)
        NUMER2 = HEIGHT(I+1) - HEIGHT(I)
        DENOM1 = HEIGHT(I) - HEIGHT(I-1)
        DENOM2 = HEIGHT(I+1) - HEIGHT(I)
        EQUAT(2) = ((NUMER1 - NUMER2) / (DENOM1 * DENOM2)) * TEMP(I)
        NUMER1 = HEIGHT(I) - HEIGHT(I-1)
        DENOM1 = HEIGHT(I+1) - HEIGHT(I-1)
        DENOM2 = HEIGHT(I+1) - HEIGHT(I)
        EQUAT(3) = (NUMER1 / (DENOM1 * DENOM2)) * TEMP(I+1)
        LPSRTE(I) = -(EQUAT(1) - EQUAT(2) + EQUAT(3))
   60 CONTINUE
C
C     THE LAST LEVEL USES THE TWO PREVIOUS LEVELS TO DETERMINE
C     THE LAPSE RATE USING THE FOLLOWING FORMULA:
C
C      LPSRTE(BOTTOM) = -((((((HEIGHT(BOTTOM)-HEIGHT(BOTTOM-1)) /
C            ((HEIGHT(BOTTOM-1)-HEIGHT(BOTTOM-2))*(HEIGHT(BOTTOM) -
C            HEIGHT(BOTTOM-2)))))) * TEMP(BOTTOM-2))
C
C             -((((HEIGHT(BOTTOM)-HEIGHT(BOTTOM-2)) /
C             ((HEIGHT(BOTTOM-1)-HEIGHT(BOTTOM-2))*(HEIGHT(BOTTOM) -
C             HEIGHT(BOTTOM-1))))) * TEMP(BOTTOM-1))
C
C             +((((2*(HEIGHT(BOTTOM)-HEIGHT(BOTTOM-1)) +
C             (HEIGHT(BOTTOM-1)-HEIGHT(BOTTOM-2))) / ((HEIGHT(BOTTOM)
C             -HEIGHT(BOTTOM-2))*(HEIGHT(BOTTOM)-HEIGHT(BOTTOM-1))))))
C              * TEMP(BOTTOM))
C
C     WHICH IS IMPLEMENTED AS FOLLOWS:
C
      NUMER1 = HEIGHT(BOTTOM) - HEIGHT(BOTTOM-1)
      DENOM1 = HEIGHT(BOTTOM-1) - HEIGHT(BOTTOM-2)
      DENOM2 = HEIGHT(BOTTOM) - HEIGHT(BOTTOM-2)
      EQUAT(1) = (NUMER1 / (DENOM1 * DENOM2)) * TEMP(BOTTOM-2)
      NUMER1 = HEIGHT(BOTTOM) - HEIGHT(BOTTOM-2)
      DENOM1 = HEIGHT(BOTTOM-1) - HEIGHT(BOTTOM-2)
      DENOM2 = HEIGHT(BOTTOM) - HEIGHT(BOTTOM-1)
      EQUAT(2) = (NUMER1 / (DENOM1 * DENOM2)) * TEMP(BOTTOM-1)
      NUMER1 = 2*(HEIGHT(BOTTOM) - HEIGHT(BOTTOM-1))
      NUMER2 = HEIGHT(BOTTOM-1) - HEIGHT(BOTTOM-2)
      DENOM1 = HEIGHT(BOTTOM) - HEIGHT(BOTTOM-2)
      DENOM2 = HEIGHT(BOTTOM) - HEIGHT(BOTTOM-1)
      EQUAT(3) = ((NUMER1 + NUMER2) / (DENOM1 * DENOM2)) * TEMP(BOTTOM)
      LPSRTE(BOTTOM) = -(EQUAT(1) - EQUAT(2) + EQUAT(3))
C
      RETURN
      END
