      FUNCTION THICKNESS(PRSTOP, PRSARY, TMPARY, DPTARY, N)
      IMPLICIT NONE
C
C**********************************************************************
C
C   NAME- THICKNESS
C
C   LANGUAGE- FORTRAN     TYPE- REAL FUNCTION  COMPUTER- AIX
C
C   VERSION- 1.0    DATE- 10/16/06   PROGRAMMER- FRANK TILLEY (IMSG)
C       ADAPTED FROM FUNCTION Z.
C
C   FUNCTION-  THIS FUNCTION RETURNS THE THICKNESS OF A LAYER
C              BOUNDED BY PRESSURE PRSARY(1) AT THE BOTTOM AND
C              PRESSURE PRSTOP AT THE TOP.
C              THE ALGORITHM INVOLVES NUMERICAL INTEGRATION
C              OF THE HYPSOMETRIC EQUATION FROM PRSARY(1) TO PRSTOP.
C              THE VALUES FOR PRSARY(X), TMPARY(X), AND
C              DPTARY(X) ARE ENTRIES FROM A SOUNDING TABLE OF N
C              ELEMENTS WITH A MAXIMUM OF 100 ELEMENTS.
C              IF N IS LARGER THAN 100 THEN THE VALUE OF THICKNESS IS
C              SET TO A MISSING VALUE (R4SPV).  IF ANY TEMPERATURE
C              VALUES BETWEEN PRSARY(I) AND PRSTOP ARE FOUND TO BE
C              MISSING, THEN R4SPV IS RETURNED FOR THE VALUE OF THICKNESS.
C              R4SPV IS ALSO RETURNED IF PRSTOP IS A MISSING VALUE.
C
C   CALLING SEQUENCE PARAMETERS-
C     VARIABLE          TYPE I/O  FUNCTION
C     --------          ---- ---  --------
C     PRSTOP            REAL I    PRESSURE (MB) AT THE TOP.
C     PRSARY            REAL I    VECTOR OF LENGTH N CONTAINING
C                                 PRESSURES (MB).
C     TMPARY            REAL I    VECTOR OF LENGTH N CONTAINING
C                                 TEMPERATURE (C) AT EACH OF THE
C                                 PRESSURE LEVELS IN THE VECTOR OF
C                                 PRESSURES.
C     DPTARY            REAL I    VECTOR OF LENGTH N CONTAINING DEW
C                                 POINT TEMPERATURES(C) AT EACH OF THE
C                                 PRESSURE LEVELS IN THE VECTOR OF
C                                 PRESSURES.
C     N                 INT  I    NUMBER OF LEVELS.
C
C
      INTEGER    N
      REAL   PRSTOP, PRSARY(N), TMPARY(N), DPTARY(N), THICKNESS
C
C   COMMON AREA PARAMETERS-
C
      INCLUDE 'RMPHSCON'
C
C   LOCAL PARAMETERS-
C     VARIABLE          TYPE    DESCRIPTION
C     --------          ----    --------
C     A1                REAL    VIRTUAL TEMPERATURE ON I+1 LEVEL.
C     A2                REAL    VIRTUAL TEMPERATURE ON I LEVEL.
C     C1                REAL    0.001*(1./EPS-1.) WHERE EPS = .62197
C                               IS THE RATIO OF THE MOLECULAR WEIGHT
C                               OF WATER TO THAT OF DRY AIR.
C                               THE FACTOR 1000. CONVERTS THE
C                               MIXING RATIO W FROM G/KG TO A
C                               DIMENSIONLESS RATIO.
C     C2                REAL    R/(2.*G) WHERE R IS THE GAS CONSTANT
C                               FOR DRY AIR (287 KG/JOULE/DEG K) AND
C                               G IS THE ACCELERATION DUE TO THE
C                               EARTH'S GRAVITY (9.8 M/S**2).  THE
C                               FACTOR OF 2 IS USED IN AVERAGING TWO
C                               VIRTUAL TEMPERATURES.
C     CMPDPT            REAL    INTERPOLATED DEW POINT TEMPERATURE
C                               AT PRSTOP - RETURNED FROM PRSINT.
C     CMPTMP            REAL    INTERPOLATED TEMPERATURE AT
C                               PRSTOP - RETURNED FROM PRSINT.
C     I                 INT     LOOP CONTROL OF LEVELS.
C     M                 INT     LOOP INDEX.
C     TI                REAL    TEMPORARY STORAGE FOR CURRENT LEVEL
C                               DEW POINT TEMPERATURE.
C     TJ                REAL    TEMPORARY STORAGE FOR NEXT LEVEL DEW
C                               POINT.
C     TK                REAL    VECTOR OF LENGTH 100 CONTAINING
C                               KELVIN TEMPERATURE EQUIVALENTS
C                               OF TEMPERATURES IN VECTOR TMPARY.
C
C   FUNCTIONS USED- DEWPOINT, PRSINT, TDEW2MR
C
      REAL       DEWPOINT, PRSINT, TDEW2MR
C
C   INTRINSIC FUNCTIONS USED- ALOG, ANINT
C
      INTEGER   I,M
      REAL      TI, TJ, C1, C2, A1, A2, TK(100), CMPTMP, CMPDPT
C
      INCLUDE 'SPVAL'
C
C   EXIT STATES-  THICKNESS RETURNS A MISSING VALUE IF THE THICKNESS CANNOT
C                 BE CALCULATED OR IF 'PRSTOP' IS A MISSING VALUE.
C
C   COMMENTS-
C   REFERENCE:    ALGORITHMS FOR THERMODYNAMIC CALCULATIONS BY
C                 THOMAS W. SCHLATTER AND DONALD V. BAKER, (1981)
C                 PROFS PROGRAM OFFICE,
C                 NOAA ENVIRONMENTAL RESEARCH LAB.
C                 BOULDER,COLORADO.
C
C   09/23/93      THE ORIGINAL VERSION OF THIS FUNCTION INCORRECTLY
C                 USED THE VIRTUAL TEMPERATURE OF THE BOUNDARY LEVEL
C                 ABOVE THE TOP PRESSURE LEVEL WHEN COMPUTING THE
C                 FINAL GEOPOTENTIAL HEIGHT.  THE FUNCTION WAS
C                 CORRECTED TO USE THE VIRTUAL TEMPERATURE AT
C                 THE TOP PRESSURE LEVEL INSTEAD.
C
C**********************************************************************
C
      LOGICAL   ISR4SPV
      PARAMETER(C1=0.0006078, C2=14.64285)
C
C     IF PRSTOP IS A MISSING VALUE SET THICKNESS TO A MISSING VALUE AND
C     RETURN.
C
      IF (ISR4SPV(PRSTOP)) THEN
         THICKNESS = R4SPV
         RETURN
      END IF
C
C     IF THE NUMBER OF LEVELS EXCEEDS 100, SET THICKNESS TO A MISSING VALUE AND
C     RETURN.
C
      IF (N.GT.100) THEN
         THICKNESS = R4SPV
         RETURN
      END IF
C
C  CONVERT TEMPERATURES FROM CELSIUS TO KELVIN.
C
      DO 10 M=1,N
        IF (.NOT. ISR4SPV(TMPARY(M))) THEN
           TK(M)=TMPARY(M)+KELVIN
        ELSE
           TK(M)=R4SPV
        ENDIF
   10 CONTINUE
C
      THICKNESS = 0.0
C
C   IF THE PRESSURE AT THE LEVEL FOR WHICH GEOPOTENTIAL
C   HEIGHT IS TO BE COMPUTED DOES NOT FALL BETWEEN THE
C   FIRST AND LAST PRESSURES IN THE PRESSURES VECTOR,
C   RETURN A MISSING VALUE FOR Z.
C
      IF ((ANINT(PRSARY(N)).GT.ANINT(PRSTOP)) .OR.
     +    (ANINT(PRSARY(1)).LT.ANINT(PRSTOP))) THEN
         THICKNESS=R4SPV
         RETURN
      ENDIF
C
C   LOOP THROUGH THE ARRAY ARGUMENTS DIMENSIONED N, ALWAYS
C   CHECKING THE NEXT LEVEL TO SEE THAT IT DOES NOT EXCEED
C   THE PRESSURE AT THE TOP.  NOTE THAT EVEN THOUGH I+1 IS
C   USED AS A SUBSCRIPT, THE HIGHEST ARRAY ELEMENT IT WILL
C   REACH IS N, AS COMMENTED BELOW.
C
      DO 20 I = 1,N
C
C  IN CASE OF A MISSING TEMPERATURE, RETURN A MISSING VALUE FOR
C  Z.
C
         IF ((ISR4SPV(TK(I))).OR.(ISR4SPV(TK(I+1)))) THEN
            THICKNESS = R4SPV
            RETURN
         ENDIF
C
C  IN CASE OF A MISSING DEW POINT, CALCULATE GEOPOTENTIAL HEIGHT BY
C  ASSUMING A RELATIVE HUMIDITY OF 10%.
C
         TI=DPTARY(I)
         TJ=DPTARY(I+1)
         IF(ISR4SPV(TI))  TI = DEWPOINT(TMPARY(I),10.)
         IF(ISR4SPV(TJ))  TJ = DEWPOINT(TMPARY(I+1),10.)
C
C  BREAK OUT OF THE LOOP WHEN THE PRESSURE AT THE LEVEL JUST ABOVE
C  THE CURRENT LEVEL IS LESS THAN THE TOP PRESSURE FOR WHICH
C  GEOPOTENTIAL HEIGHT IS TO BE CALCULATED.  AT THIS POINT
C  THE VALUE OF Z REFLECTS THE GEOPOTENTIAL HEIGHT OF THE CURRENT
C  LEVEL.
C
         IF (ANINT(PRSTOP).GE.ANINT(PRSARY(I+1))) GO TO 30
         A1 = TK(I+1)*(1. +C1*TDEW2MR(PRSARY(I+1),TJ))
         A2 = TK(I)*(1.+C1*TDEW2MR(PRSARY(I),TI))
         THICKNESS = THICKNESS + C2*(A1+A2)*(ALOG(PRSARY(I)/PRSARY(I+1)))
C
   20 CONTINUE
   30 CONTINUE
C
C  USE FUNCTION PRSINT TO COMPUTE AN INTERPOLATED TEMPERATURE
C  AND DEW POINT TEMPERATURE (C) FOR THE PRESSURE AT THE TOP.
C
      CMPTMP=PRSINT(TK(I+1),TK(I),PRSARY(I+1),PRSARY(I),PRSTOP)
      CMPDPT=PRSINT(TJ,TI,PRSARY(I+1),PRSARY(I),PRSTOP)
C
C  IF AN INTERPOLATED TEMPERATURE CANNOT BE COMPUTED, RETURN
C  Z AS A MISSING VALUE.
C
      IF (ISR4SPV(CMPTMP)) THEN
         THICKNESS = R4SPV
         RETURN
      ENDIF
C
C  COMPUTE THE VIRTUAL TEMPERATURES AT THE TOP PRESSURE LEVEL
C  AND AT THE NEXT LOWER BOUNDARY LEVEL.
C
      A1 = CMPTMP*(1.+C1*TDEW2MR(PRSTOP,CMPDPT))
      A2 = TK(I)*(1.+C1*TDEW2MR(PRSARY(I),TI))
C
C  USE THE VIRTUAL TEMPERATURES AT THE TOP AND AT THE NEXT LOWER
C  BOUNDARY LAYER TO COMPUTE THE GEOPOTENTIAL HEIGHT DELTA BETWEEN
C  THE BOUNDARY LAYER AND THE TOP PRESSURE LEVEL. ADD THIS VALUE TO
C  THE ACCUMULATED GEOPOTENTIAL HEIGHT OF THE LOWER BOUNDARY LAYER.
C
      THICKNESS = THICKNESS+C2*(A1+A2)*(ALOG(PRSARY(I)/PRSTOP))
C
      RETURN
C
      END
