C$$$   SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    HB2P        INTERPOLATE FROM HYB-B TO ISOBARIC GRIDS
C
C   PRGMMR:  BENJAMIN, STAN ORG: ERL/PROFS      DATE: 93-01-12
C
C ABSTRACT:  INTERPOLATE GRIDS FROM HYBRID-B LEVELS (NO INTERSECTING
C            SURFACES) TO ISOBARIC LEVELS
C
C PROGRAM HISTORY LOG:
C    93-01-12      S. BENJAMIN     ORIGINAL VERSION
C
C USAGE:           CALL   HB2P (g3,XT,
C     1            H,T,RHP,UP,VP,XP,NX3D,PLEV,NZP)
C
C   INPUT ARGUMENT LIST:
C    g3        - real
C       contains the following
C     PT       - REAL     HYBRID-B PRESSURE (Pa)
C     TH       - REAL     HYBRID-B VIRTUAL POTENTIAL TEMP (K)
C     M        - REAL     HYBRID-B height (m)
C     UT       - REAL     HYBRID-B U WIND COMPONENT (M/S)
C     VT       - REAL     HYBRID-B V WIND COMPONENT (M/S)
C     vertical velocity, cloud variables, TKE

C     XT       - REAL     other HYBRID-B variables
C     RHB      - REAL     HYBRID-B Relative humidity (0-1)
C     NX3D     - INTEGER  NUMBER OF other hybrid-b variables
C     PLEV     - REAL     ISOBARIC LEVELS IN MB
C     NZP      - INTEGER  NUMBER OF LEVELS TO INTERPOLATE TO
C   OUTPUT ARGUMENT LIST:
C     H        - REAL     ISOBARIC HEIGHT (M)
C     T        - REAL     ISOBARIC TEMPERATURE (K) (NON-VIRTUAL)
C     RHP      - REAL     ISOBARIC RELATIVE HUMIDITY (0-1)
C     UP       - REAL     ISOBARIC U WIND COMPONENT (M/S)
C     VP       - REAL     ISOBARIC V WIND COMPONENT (M/S)
C     XP       - REAL     other ISOBARIC variables
C
C REMARKS:  SEE BENJAMIN ET AL., 1993, AMS CONFERENCE ON AVIATION
C    WEATHER FORECASTING, FOR DISCUSSION OF HYBRID-B COORDINATES.
C    THE MAXIMUM NUMBER OF ISOBARIC LEVELS IS SET AT 37.  THIS CAN
C    BE INCREASED, IF NECESSARY.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN-77 + EXTENSIONS
C   MACHINE: DEC - VAX, VMS
C
C$$$
        SUBROUTINE HB2P (g3,XT,RHB,
     1            H,T,RHP,UP,VP,XP,NX3D,PLEV,NZP)

C ***   This program converts a hybrid analysis/forecast in hybrid-b
C       coordinates (non-intersecting levels) into
C       typical isobaric variables on isobaric surfaces.


C ***   NZP is the number of pressure levels to which to interpolate.  
C       The maximum value
C       for NZP is 39, which is the number of levels
C       spaced 25 mb apart between 1050 and 100 mb.  If you need
C       to interpolate to more than 39 isobaric levels, change
C       NZP25_P in MAPSCON and recompile.

C       NZTN_P is the number of hybrid levels.  An extra level
C       is needed for a higher isentropic level (500 K) with which 
C       to extrapolate to higher isobaric levels.

C       NX3D is the number of extra 3d variables to be interpolated
C       to isobaric surfaces.  Currently, these can be used for cloud
C       and TKE variables passed in thru g3.

        IMPLICIT NONE
c       INCLUDE 'MAPSCON'

        INTEGER
     1   I,J,K,NZP,N,NX3D,IVAR
        integer nx_p,ny_p,nztn_p,nvar3_p
        parameter (nx_p=151,ny_p=113,nztn_p=40,nvar3_p=14)
        integer p_p,theta_p,u_p,v_p,h_p
        parameter (p_p=1,theta_p=3,u_p=5,v_p=6,h_p=2)
        integer nzp25_p,nzp_p
        parameter (nzp25_p=37,nzp_p=8)

C *** HYBRID VARIABLES
        REAL     G3 (NX_P,NY_P,NZTN_P,NVAR3_P)
C pressure - Pa
C virtual pot temp - K
C height - m
C u component of wind - m/s
C v component of wind - m/s
C other 3-d variables (e.g., pressure vertical velocity)
     1  ,RHB(NX_P,NY_P,NZTN_P)
C relative humidity (0-1) on hybrid levels
     1  ,xt(NX_P,NY_P,NZTN_P)
C extra variable - potential vorticity

       REAL
     1   DUM  (NX_P,NY_P,NZTN_P)
     1  ,DUMMY (NX_P,NY_P)
     1  ,TEMP,Q,PM1
     1  ,EX1
C R*gamma/Cp
     1  ,GAMMA
     1  ,GAM,GAMD
     1  ,PSFC,TSFC,ZSFC
     1  ,T1,T5,Z1,Z5
C sfc values
     1  ,C (4,NZTN_P+1)
C work array for VERINT
     1  ,XI(NZTN_P+1)
C work array for VERINT

C *** THETA VARIABLES
        REAL
     1   EXNER (NX_P,NY_P,NZTN_P)

C *** PRESSURE VARIABLES
        REAL     H     (NX_P,NY_P,NZP25_P)
C height - m
     1  ,T     (NX_P,NY_P,NZP25_P)
C temperature (non-virtual)
C                                                 - K
     1  ,RHP   (NX_P,NY_P,NZP25_P)
C relative humidity
C                                                 - 0-1 scale
     1  ,UP    (NX_P,NY_P,NZP25_P)
C u component of wind - m/s
     1  ,VP    (NX_P,NY_P,NZP25_P)
C v component of wind - m/s
     1  ,XP (NX_P,NY_P,NZP25_P,NX3D)
C other 3-d variables, e.g., vertical velocity

C *** PLEV are the pressure levels in mb 
        REAL     PLEV  (NZP)
     1  ,EXLEV (NZP25_P)

C *** SMOOTHING VARIABLES
        REAL PMLEV (NZP_P)
        DATA PMLEV / 85000., 70000., 50000., 40000.,
     1               30000., 25000., 20000., 15000.  /

        INTEGER  KSM   (NZP25_P)
C         - closest mandatory level
        INTEGER  NSM (5, NZP_P)
C         - look-up table of recommended number of smoothing
C           passes, function of level and variable

C                    H    T   RH   U   V
        DATA NSM /   3,   2,   1,  3,  3,
     1               2,   1,   1,  3,  3,
     1               2,   1,   1,  2,  2,
     1               2,   1,   1,  1,  1,
     1               2,   1,   1,  1,  1,
     1               2,   1,   1,  1,  1,
     1               2,   1,   1,  1,  1,
     1               3,   1,   1,  2,  2   /

C        DATA     EX1     /0.1903643/ 
C R*gamma/g
        DATA     GAMMA   /0.0065/
C 6.5 K/km
        DATA     GAMD    /0.0100/
C 10  K/km

        REAL    CPD_P, ROVCP_P, G0_P, RD_P
        PARAMETER (
     + CPD_P   = 1004.686,
C J/(kg*K)     Specific heat at constant pressure

C                <CPD> = 3.5*<RD>
     + ROVCP_P = .285714,
C nd           <RD>/<CPD>
     +  G0_P    = 9.80665
C m/(s*s)   Sea-level acceleration of gravity
     +  ,RD_P = 287.04
     1  )

        INTEGER IG,K5
C IG - LEVEL USED FOR EXTRAPOLATING BELOW SURFACE

C        TYPE *,' ENTER LEVEL FOR EXTRAPOLATING BELOW SURFACE'
C        ACCEPT *,IG
        IG = 4

C        TYPE *,' ENTER TOP LEVEL FOR COMPUTING LAPSE RATE'
C        ACCEPT *,K5
        K5 = 5

         DO 25 K = 1,NZTN_P
         DO 25 J = 1,NY_P
         DO 25 I = 1,NX_P
          EXNER(I,J,K) = CPD_P*(g3(I,J,K,p_p)/100000.)**ROVCP_P
25      CONTINUE

        do 27 k=1,nztn_p
          print *,' exner,pt,th,rhb',exner(15,15,k),g3(15,15,k,p_p),
     1      g3(15,15,k,theta_p),rhb(15,15,k)
27      continue

        DO  10 K = 1,NZP
          EXLEV(K) = CPD_P*(PLEV(K)/100000.)**ROVCP_P
          PRINT *,' K,P,EX',K,PLEV(K),EXLEV(K)
10      CONTINUE

        CALL VERINT (0,1,EXNER,g3(1,1,1,u_p),DUM,NX_P,NY_P,NZTN_P,C,XI,
     1  UP,NZP,EXLEV)

        CALL VERINT (0,1,EXNER,g3(1,1,1,v_p),DUM,NX_P,NY_P,NZTN_P,C,XI,
     1  VP,NZP,EXLEV)

        DO 11 IVAR = 1,NX3D-1
          CALL VERINT (0,1,EXNER,g3(1,1,1,6+IVAR),DUM,NX_P,NY_P,
     1      NZTN_P,C,XI,XP(1,1,1,IVAR),NZP,EXLEV)
11      CONTINUE

        CALL VERINT (0,1,EXNER,RHB,DUM,NX_P,NY_P,NZTN_P,C,XI,
     1  RHP,NZP,EXLEV)

        CALL VERINT (0,1,EXNER,xt,DUM,NX_P,NY_P,NZTN_P,C,XI,
     1  xp(1,1,1,nx3d),NZP,EXLEV)

        CALL VERINT (0,1,EXNER,g3(1,1,1,theta_p),DUM,NX_P,NY_P,NZTN_P,
     1    C,XI,T,NZP,EXLEV)

        CALL VERINT (0,1,EXNER,g3(1,1,1,h_p),DUM,NX_P,NY_P,NZTN_P,
     1    C,XI,H,NZP,EXLEV)

        PRINT *,' H(15,15,5)',H(15,15,5)
        PRINT *,' T(15,15,3)',T(15,15,3)
        PRINT *,' RHP(15,15,5)',RHP(15,15,5)
        PRINT *,' UP(15,15,5)',UP(15,15,5)
        PRINT *,' VP(15,15,5)',VP(15,15,5)
        DO 12 IVAR = 1,NX3D
        PRINT *,' XP(15,15,5),IVAR',XP(15,15,5,IVAR)
12      CONTINUE

C________________________________________________________
C________________________________________________________
C     Extrapolate and interpolate near ground
C________________________________________________________
C________________________________________________________
           DO 50 J = 1,NY_P
           DO 50 I = 1,NX_P
            PSFC = 100000.*(EXNER(I,J,IG)/CPD_P)**3.5
            TSFC = EXNER(I,J,IG)*g3(I,J,IG,theta_p)/CPD_P
            ZSFC = g3(I,J,IG,h_p)

Cxxxxxx ERROR, should be theta at level 1
            T1 = EXNER(I,J,1)*g3(I,J,IG,theta_p)/CPD_P
            Z1 = g3(I,J,1,h_p)

            T5 = EXNER(I,J,K5)*g3(I,J,K5,theta_p)/CPD_P
            Z5 = g3(I,J,K5,h_p)

            GAM = (T1-T5)/(Z5-Z1)
            GAM = MIN(GAMD,MAX(GAM,GAMMA))
            EX1 = RD_P*GAM/G0_P
         DO 50 K = 1,NZP

C---------------------------------
C  The interpolation is considered in 3 layers: A,B,C
C    A = above IGth hybrid-b level (about 15-25 hPa above ground)
C    B = between IGth level and model surface
C    C = below model surface (beneath ground)
C
C         ht      temp    RH     u    v    x
C       ------   ------  ----- ----- ---- -----
C  A    interp   interp  inter inter inter inter
C - - - - - - - - - - - - - - - - - - - - - - - -
C  B    reduce   reduce  inter inter inter inter
C _______________________________________________
C / / / / / / / / / surface / / / / / / / / / / /
C  C    reduce   reduce  =k=1  =k=1  =k=1  =k=1
C
C - layer A
            IF (EXLEV(K).LT.EXNER(I,J,IG)) THEN

C --- get T from theta

            T(I,J,K) = T(I,J,K)*EXLEV(K)/CPD_P
            TEMP         = T(I,J,K)


C - layer B
          ELSE IF (EXLEV(K).GT.EXNER(I,J,IG) .AND.
     1             EXLEV(K).LT.EXNER(I,J,1)) THEN
C --- we need to do some extrapolating downward at this grid point.
C       Set RH, u, v = sfc values for isobaric levels beneath ground.
C in pascals
            T(I,J,K) = TSFC*(PLEV(K)/PSFC)**EX1
            H(I,J,K) = ZSFC - (T(I,J,K)-TSFC)/GAM
          ELSE
C - layer C   --  IF (EXLEV(K).GT.EXNER(I,J,1)) THEN
            T(I,J,K) = TSFC*(PLEV(K)/PSFC)**EX1
            H(I,J,K) = ZSFC - (T(I,J,K)-TSFC)/GAM
            UP(I,J,K) = g3(I,J,1,u_p)
            VP(I,J,K) = g3(I,J,1,v_p)
            RHP(I,J,K)=RHB(I,J,1)
            DO 49 IVAR=1,NX3D-1
            XP(I,J,K,IVAR) = g3(I,J,1,6+IVAR)
49          CONTINUE
            xp(i,j,k,nx3d) = xt(i,j,1)
          END IF
c      if (k.eq.1.and.(t(i,j,k).lt.200..or.t(i,j,k).gt.320.))
c    1    print *,' bad reduction',i,j,k,t(i,j,k),h(i,j,k),
c    1     tsfc,zsfc,psfc,gam,gamd,gamma
50      CONTINUE

C *** De-virtualize virtual temperatures 
         DO 55 K = 1,NZP
         DO 55 J = 1,NY_P
         DO 55 I = 1,NX_P
          CALL TV2TQ (T(I,J,K),RHP(I,J,K)/100.,PLEV(K)/100.,
     1  TEMP,Q)
c      if (k.eq.1.and.(temp.lt.200..or.temp.gt.320.)) then
c         print *,' bad tv2tq',i,j,k,temp,t(i,j,k),rhp(i,j,k)
c         do k1=1,25
c          print *,t(i,j,k1),rhp(i,j,k1),pt(i,j,k1),th(i,j,k1)
c         end do
c       end if
          T(I,J,K) = TEMP
55      CONTINUE

C --- Find closest mandatory levels to each isobaric level
         DO 58 K = 1,NZP
           CALL FNDLVL (PLEV(K),PMLEV,NZP_P,KSM(K))
58       CONTINUE

C________________________________________________________
C________________________________________________________
C  Extrapolate to upper isobaric levels, if necessary
C    Check from 175 to 100 mb.
C________________________________________________________
C________________________________________________________

        do 52 k=nzp-3,nzp
        do 52 j=1,ny_p
        do 52 i=1,nx_p
          if (t(i,j,k).lt.0.) then
            t(i,j,k) = g3(i,j,nztn_p,theta_p)*exner(i,j,nztn_p)/cpd_p
            h(i,j,k) = g3(i,j,nztn_p,h_p) + rd_p*t(i,j,k)/g0_p
     +           *alog(g3(i,j,nztn_p,p_p)/plev(k))
            rhp(i,j,k) = rhb(i,j,nztn_p)
            up (i,j,k) = g3 (i,j,nztn_p,u_p)
            vp (i,j,k) = g3 (i,j,nztn_p,v_p)
            do 515 ivar = 1,nx3d-1
              xp(i,j,k,ivar) = g3(i,j,nztn_p,6+ivar)
515         continue
            xp(i,j,k,nx3d) = xt(i,j,nztn_p)
          end if
52      continue


C --- Smooth by number of passes for recommended level
C       at nearest mandatory level (NSM)

         DO 60 K = 1,NZP
          DO 591 N=1,NSM(1,KSM(K))
            CALL SMOOTH(H(1,1,K),DUMMY,NX_P,NY_P,0.5)
591       CONTINUE
          DO 592 N=1,NSM(2,KSM(K))
            CALL SMOOTH(T(1,1,K),DUMMY,NX_P,NY_P,0.5)
592       CONTINUE
          DO 593 N=1,NSM(3,KSM(K))
            CALL SMOOTH(RHP(1,1,K),DUMMY,NX_P,NY_P,0.5)
            DO 5925 IVAR=1,NX3D
            CALL SMOOTH(XP(1,1,K,IVAR),DUMMY,NX_P,NY_P,0.5)
5925        CONTINUE
593       CONTINUE
          DO 594 N=1,NSM(4,KSM(K))
            CALL SMOOTH(UP(1,1,K),DUMMY,NX_P,NY_P,0.5)
            CALL SMOOTH(VP(1,1,K),DUMMY,NX_P,NY_P,0.5)
594       CONTINUE

60      CONTINUE

         DO 70 K = 1,NZP
         DO 70 J = 1,NY_P
         DO 70 I = 1,NX_P
            RHP(I,J,K) = MAX(0.,MIN(100.,RHP(I,J,K)))
70      CONTINUE

        RETURN
        END

C$$$  SUBPROGRAM DOCUMENTATION BLOCK 
C                .      .    .                                       . 
C SUBPROGRAM:    VERINT      PERFORM VERTICAL INTERPOLATION
C   PRGMMR: STAN BENJAMIN    ORG: FSL/PROFS  DATE: 90-09-19 
C 
C ABSTRACT: COMPUTE VERTICAL HERMITE OR LINEAR INTERPOLATION. 
C 
C PROGRAM HISTORY LOG: 
C   
C 
C USAGE:    CALL VERINT (ITYPE, CALCDF, SUPORT, F, DF, 
C                        IDIN, JDIN, KDIN1, C, XI, 
C                        FNEW,  KDIN2, NEWLEV)
C   INPUT ARGUMENT LIST: 
C--MODEV
C     ITYPE    - INTEGER    TYPE OF INTERPOLATION SCHEME.   
C                           0 = LINEAR, 1 = HERMITE
C     CALCDF   - INTEGER    DECIDE WHETHER USE THE DERIVATIVE VALUE 
C                FROM OUTSIDE OF ROUTINE OR NOT.
C                0 = USE,  1 = CALCULATE IN THE ROUTINE. 
C--QI
C     SUPORT   - REAL ARRAY   SUPORT(IDIN,JDIN,KDIN1)
C                SUPPORT POINTS IN NEW VERTICAL 
C                COORDINATES ON OLD LEVEL.  E.G. THETA
C--FI
C     F        - REAL ARRAY   F(IDIN,JDIN,KDIN1)
C                VALUES AT SUPPORT POINT TO INTERPOLATE.
C--DFI
C     DF       - REAL ARRAY   DF(IDIN,JDIN,KDIN1)
C                DERIVATIVE OF F FOR HERMITE INTERPOLATION.  
C     IDIN     - INTEGER   NUMBER OF GRID POINTS IN X-DIRECTION.
C     JDIN     - INTEGER   NUMBER OF GRID POINTS IN Y-DIRECTION.
C     KDIN1    - INTEGER   NUMBER OF GRID POINTS IN VERTICAL LEVEL.
C     C        - REAL ARRAY   C(4,KDIN1)
C                COEFFICIENTS OF CUBIC POLYNOMIAL BETWEEN 2 POINTS.
C     XI       - REAL ARRAY   XI(KDIN1)
C                TEMPORARY STORAGE VALUE TO BE INTERPOLATED.
C--KDOU
C     KDIN2    - INTEGER    NUMBER OF GRID POINTS IN VERTICAL LEVEL
C                ON NEW COORDINATES(NEWLEV).
C--QKO
C     NEWLEV   - REAL ARRAY   NEWLEV(KDIN2)
C                NEW COORDINATES LEVELS TO EVALUATE
C                INTERPOLATING POLYNOMIALS AT.
C 
C   OUTPUT ARGUMENT LIST:   
C--FO
C     FNEW     - REAL ARRAY   FNEW(IDIN,JDIN,KDIN2)
C                INTERPOLATED VALUES OF F ON NEW COORDINATES LEVEL.
C 
C REMARKS: THIS PROGRAM WORKS FOR ONLY MONOTONICALLY DECREASING
C   SUPPORT POINT.  ALSO, WHEN THE DERIVATIVE VALUE OF F (=DF) ARE
C   USED FROM PASSING ARGUEMENT ( WHEN CALCDF = 0), THE VALUES OF
C   DF CORRESPOND TO THE MONOTONICALLY DECREASING SUPPORT POINTS.
C 
C ATTRIBUTES: 
C   LANGUAGE: FORTRAN-77
C   MACHINE:  VAX-8800
C$$$ 
C
      SUBROUTINE VERINT (ITYPE, CALCDF, SUPORT, F, DF, 
     -                   IDIN, JDIN, KDIN1, C, XI, 
     -                   FNEW,  KDIN2, NEWLEV)
C                
      IMPLICIT NONE
      INTEGER    ITYPE, CALCDF, IDIN, JDIN, KDIN1, KDIN2
      REAL       SUPORT, F, FNEW, NEWLEV
      REAL       DF, XI, C

      DIMENSION  SUPORT(IDIN,JDIN,KDIN1),  F(IDIN,JDIN,KDIN1),
     -           DF(IDIN,JDIN,KDIN1)
      DIMENSION  XI(KDIN1),    C(4,KDIN1)
      DIMENSION  FNEW(IDIN,JDIN,KDIN2),
     -           NEWLEV(KDIN2)
      INTEGER    LIN, HERM
      PARAMETER  ( LIN = 0, HERM = 1 )
C
C
      IF (ITYPE .EQ. LIN) THEN
        CALL LININT ( SUPORT, CALCDF, F, DF, IDIN, JDIN, KDIN1, 
     -                FNEW,  KDIN2, NEWLEV, XI)

      ELSEIF (ITYPE .EQ. HERM) THEN
        CALL HERINT ( SUPORT, CALCDF, F, DF, IDIN, JDIN, KDIN1, 
     -                FNEW,  KDIN2, NEWLEV, C, XI)

      ELSE
        PRINT*, ' INVALID INTERPOLATION TYPE. ITYPE =', ITYPE
      ENDIF

      RETURN
       END


C$$$   SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    LININT      PERFORM LINEAR INTERPOLATION.
C   PRGMMR:  BENJAMIN, STAN ORG: ERL/PROFS      DATE: 93-01-18
C
C ABSTRACT:  PERFORM LINEAR INTERPOLATION.
C
C PROGRAM HISTORY LOG:
C    88/05/31       S. BENJAMIN     ORIGINAL VERSION
C
C USAGE:   CALL LININT ( SUPORT, CALCDF, F, DF, IDIN, JDIN, KDIN1,
C                    FNEW,  KDIN2, NEWLEV, XI)
C
C   INPUT ARGUMENT LIST:
C
C   OUTPUT ARGUMENT LIST:
C
C REMARKS: NONE
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN-77 + EXTENSIONS
C   MACHINE: DEC - VAX, VMS
C
C$$$
C-------------------------------------------------------------------
C                   PERFORM LINEAR INTERPOLATION.
C-------------------------------------------------------------------
      SUBROUTINE LININT ( SUPORT, CALCDF, F, DF, IDIN, JDIN, KDIN1, 
     -                    FNEW,  KDIN2, NEWLEV, XI)
C
      IMPLICIT NONE
      INTEGER  CALCDF, IDIN, JDIN, KDIN1, KDIN2, I, J, K, KK, PFIND
      REAL     XBAR
      REAL     SUPORT(IDIN,JDIN,KDIN1),  F(IDIN,JDIN,KDIN1)
      REAL     DF(IDIN,JDIN,KDIN1),     
     -         XI(KDIN1)
C
      REAL     FNEW(IDIN,JDIN,KDIN2),
     -         NEWLEV(KDIN2) 
C
      IF (CALCDF.EQ.1)
     .THEN
C
C  CALCULATE DERIVATIVE OF F.
C
         DO 100 J=1, JDIN
         DO 100 I=1, IDIN
         DO 100 K=1, KDIN1-1
          DF(I,J,K) = (F(I,J,K+1) - F(I,J,K))/
     -                 (SUPORT(I,J,K+1) - SUPORT(I,J,K))      
  100   CONTINUE
C
         DO 110 J=1, JDIN
         DO 110 I=1, IDIN
          DF(I,J,KDIN1) = DF(I,J,KDIN1-1)
  110   CONTINUE
      ENDIF
C
C GATHER SUPPORT POINT
C 
       DO 160 J=1, JDIN
       DO 160 I=1, IDIN
         DO 140 K=1, KDIN1
          XI(K)= SUPORT(I,J,K)
  140   CONTINUE  
C
C EVALUATE INTERPOLATION.
C                                             
         DO 150 K=1, KDIN2
          XBAR = NEWLEV(K)
          KK =  PFIND(XBAR, XI, KDIN1-1)
          FNEW(I,J,K) = F(I,J,KK) + DF(I,J,KK)*(XBAR - XI(KK))
  150   CONTINUE
  160 CONTINUE
C          
      RETURN
       END


C$$$   SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    PFIND       RETURN LINEAR INTERP INDEX
C   PRGMMR:  BENJAMIN, STAN ORG: ERL/PROFS      DATE: 93-01-18
C
C ABSTRACT:  RETURNS THE INDEX VALUE AT XBAR OF THE LINEAR INTERPOLATION
C     ON N INTERVALS WITH BREAKPOINT SEQUENCE XI.
C
C PROGRAM HISTORY LOG:
C    88/05/31       S. BENJAMIN     ORIGINAL VERSION
C
C USAGE:   CALL PFIND(XBAR, XI, N)
C
C   INPUT ARGUMENT LIST:
C
C   OUTPUT ARGUMENT LIST:
C
C REMARKS: NONE
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN-77 + EXTENSIONS
C   MACHINE: DEC - VAX, VMS
C
C$$$
C-------------------------------------------------------------------
C RETURNS THE INDEX VALUE AT XBAR OF THE LINEAR INTERPOLATION 
C ON N INTERVALS WITH BREAKPOINT SEQUENCE XI.
C-------------------------------------------------------------------
      INTEGER FUNCTION PFIND (XBAR, XI, N)

      IMPLICIT NONE
      INTEGER N, I, J
      REAL    XBAR, XI(N+1)
      DATA    I /1/

      IF (XBAR .LE. XI(I)) THEN
         DO 555 J=I, N
          IF (XBAR .GT. XI(J+1))  GO TO 777
 555    CONTINUE
        J = N
      ELSE
         DO 666 J=I-1, 1, -1 
          IF (XBAR .GE. XI(J))    GO TO 777
 666    CONTINUE
        J = 1
       ENDIF

 777  PFIND = J

      RETURN
       END


C$$$   SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    HERINT      PERFORM HERIMITE INTERPOLATION
C   PRGMMR:  BENJAMIN, STAN ORG: ERL/PROFS      DATE: 93-01-18
C
C ABSTRACT:  PERFORM HERMITE INTERPOLATION
C
C PROGRAM HISTORY LOG:
C    88/05/31       S. BENJAMIN     ORIGINAL VERSION
C
C USAGE:   CALL HERINT ( SUPORT, CALCDF, F, DF,
C    1                    IDIN, JDIN, KDIN1, 
C    2                    FNEW,  KDIN2, NEWLEV, C, XI)
C
C   INPUT ARGUMENT LIST:
C
C   OUTPUT ARGUMENT LIST:
C
C REMARKS: NONE
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN-77 + EXTENSIONS
C   MACHINE: DEC - VAX, VMS
C
C$$$
C--------------------------------------------------------------------
C                     PERFORM HERMITE INTERPOLATION.
C--------------------------------------------------------------------
      SUBROUTINE HERINT ( SUPORT, CALCDF, F, DF, 
     -                    IDIN, JDIN, KDIN1, 
     -                    FNEW,  KDIN2, NEWLEV, C, XI)
C
      IMPLICIT NONE
      REAL     XBAR, PCUBIC
      INTEGER  CALCDF, IDIN, JDIN, KDIN1, KDIN2, I, J, K
      REAL     SUPORT(IDIN,JDIN,KDIN1),  F(IDIN,JDIN,KDIN1)
      REAL     DF(IDIN,JDIN,KDIN1),      
     -         C(4,KDIN1), XI(KDIN1)
C
      REAL     FNEW(IDIN,JDIN,KDIN2),
     -         NEWLEV(KDIN2) 
C
      IF (CALCDF.EQ.1) THEN
C
C CALCULATE THE DERIVATIVE OF F.
C
         DO 200 K=2, KDIN1 - 1
         DO 200 J=1, JDIN
         DO 200 I=1, IDIN
           DF(I,J,K) = (F(I,J,K+1) - F(I,J,K-1))/
     -                 (SUPORT(I,J,K+1) - SUPORT(I,J,K-1))
 200    CONTINUE

         DO 210 J=1, JDIN
         DO 210 I=1, IDIN
           DF(I,J,1) = (F(I,J,2) - F(I,J,1))/
     -                 (SUPORT(I,J,2) - SUPORT(I,J,1))
           DF(I,J,KDIN1) = (F(I,J,KDIN1) - F(I,J,KDIN1-1))/
     -                   (SUPORT(I,J,KDIN1) - SUPORT(I,J,KDIN1-1))
 210    CONTINUE
      ENDIF
C
C  GATHER SUPPORT POINT AND CALCULATE THE CUBIC COEFFICIENT
C  OF EACH INTERVAL OF F.
C
       DO 250 J=1, JDIN
       DO 250 I=1, IDIN
         DO 235 K=1,KDIN1
            XI(K) = SUPORT(I,J,K)
            C(1,K) = F(I,J,K)
            C(2,K) = DF(I,J,K)
 235    CONTINUE
        CALL CALCF( XI,C,KDIN1-1 )
         DO 240 K=1,KDIN2
          XBAR = NEWLEV(K)
          FNEW(I,J,K) = PCUBIC(XBAR,XI,C,KDIN1-1)
 240    CONTINUE
 250  CONTINUE
      RETURN
       END



C$$$   SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    CALCF       ???
C   PRGMMR:  BENJAMIN, STAN ORG: ERL/PROFS      DATE: 93-01-18
C
C ABSTRACT:  ???
C
C PROGRAM HISTORY LOG:
C    88/05/31       S. BENJAMIN     ORIGINAL VERSION
C
C USAGE:   CALL CALCF (XI, C, N)
C
C   INPUT ARGUMENT LIST:
C
C   OUTPUT ARGUMENT LIST:
C
C REMARKS:
C INPUT;
C XI(1), ..., X(I+1) STRICTLY INCREASING SEQUENCE OF BREAKPOINTS.
C C(1,I), C(2,I), VALUES AND FIRST DERIVATIVE AT XI(I), I=1,...,N+1
C
C OUTPUT;
C POLYNOMIAL COEFFICIENTS FOR C(1,I), C(2,I), C(3,I), C(4,I) OF
C INTERVAL (XI(I), XI(I+1)).  I= 1,..,N
C-------------------------------------------------------------------
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN-77 + EXTENSIONS
C   MACHINE: DEC - VAX, VMS
C
C$$$
C-------------------------------------------------------------------
C INPUT;
C XI(1), ..., X(I+1) STRICTLY INCREASING SEQUENCE OF BREAKPOINTS.
C C(1,I), C(2,I), VALUES AND FIRST DERIVATIVE AT XI(I), I=1,...,N+1
C OUTPUT;
C POLYNOMIAL COEFFICIENTS FOR C(1,I), C(2,I), C(3,I), C(4,I) OF
C INTERVAL (XI(I), XI(I+1)).  I= 1,..,N
C-------------------------------------------------------------------
      SUBROUTINE CALCF (XI, C, N)
      INTEGER N, I
      REAL    C(4,N+1), XI(N+1), DIVDF1, DIVDF3, DX

       DO 12 I=1, N
        DX = XI(I+1) - XI(I)
        DIVDF1 = (C(1,I+1) - C(1,I))/DX
        DIVDF3 = C(2,I) + C(2,I+1) - 2.0*DIVDF1
        C(3,I) = (DIVDF1 - C(2,I) - DIVDF3)/DX
        C(4,I) = DIVDF3/(DX*DX)
 12   CONTINUE
      RETURN
       END



C$$$   SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    PCUBIC      PIECEWISE CUBIC FUNCTIONS
C   PRGMMR:  BENJAMIN, STAN ORG: ERL/PROFS      DATE: 93-01-18
C
C ABSTRACT: PIECEWISE CUBIC FUNCTION 
C
C PROGRAM HISTORY LOG:
C    88/05/31       S. BENJAMIN     ORIGINAL VERSION
C
C USAGE:   PCUBIC (XBAR, XI, C, N)
C
C   INPUT ARGUMENT LIST:
C
C   OUTPUT ARGUMENT LIST:
C
C REMARKS: NONE
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN-77 + EXTENSIONS
C   MACHINE: DEC - VAX, VMS
C
C$$$
C-------------------------------------------------------------------
C RETURNS THE VALUE AT XBAR OF THE PIECEWISE CUBIC FUNCTION
C ON N INTERVALS WITH BREAKPOINT SEQUENCE XI AND COEFFICIENTS C.
C-------------------------------------------------------------------
      REAL FUNCTION PCUBIC (XBAR, XI, C, N)

      IMPLICIT NONE
      INTEGER N, I, J
      REAL  C(4, N), XBAR, XI(N+1), DX
      DATA  I /1/

      IF (XBAR .LE. XI(I)) THEN
         DO 222 J=I, N
          IF (XBAR .GT. XI(J+1))  GO TO 444
 222    CONTINUE
        J = N
      ELSE
         DO 333 J=I-1, 1, -1 
          IF (XBAR .LE. XI(J))    GO TO 444
 333    CONTINUE
        J = 1
       ENDIF

 444  I = J

      DX = XBAR - XI(I)
      PCUBIC = C(1,I) + DX*(C(2,I) + DX*(C(3,I) + DX*C(4,I)))

      RETURN
       END



C$$$   SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    TV2TQ       DEVIRTUALIZE TEMPERATURE
C
C   PRGMMR:  BENJAMIN, STAN ORG: ERL/PROFS      DATE: 93-01-12
C
C ABSTRACT:  CALCULATES TEMPERATURE AND SPECIFIC HUMIDITY FROM 
C            VIRTUAL TEMPERATURE, PRESSURE, AND RELATIVE HUMIDITY
C
C PROGRAM HISTORY LOG:
C   92       S. BENJAMIN     ORIGINAL VERSION
C   94-05-16 S. BENJAMIN     IMPROVE ACCURACY BY USING SPECIFIC 
C                             HUMIDITY INSTEAD OF MIXING RATIO.
C                             ALSO HAS FASTER CONVERGENCE.
C                             ADD ORDER OF ACCURACY TO 0.6078
C                             CONSTANT.
C   94-05-23 S. BENJAMIN     LIMIT ALGORITHM TO 3 ITERATIONS
C                             FOR MORE SPEED.
C   95-01-02 S. BENJAMIN     NEW VERSION
C   95-07-28 S. Benjamin     Limit number of iterations
C
C USAGE:           CALL TV2TQ(TV, RH, P, T, Q)
C
C   INPUT ARGUMENT LIST:
C   TV          - REAL           VIRTUAL TEMPERATURE IN KELVIN
C   RH          - REAL           RELATIVE HUMIDITY
C   P           - REAL           PRESSURE IN MILIBARS
C
C   OUTPUT ARGUMENT LIST:
C   T           - REAL           TEMPERATURE IN KELVIN
C   Q           - REAL           SPECIFIC HUMIDITY (G/G)
C
C REMARKS: NONE
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN-77 + EXTENSIONS
C   MACHINE: DEC - VAX, VMS
C

C$$$
        SUBROUTINE TV2TQ(TV,RH,P,T,Q)
C
C   THIS ROUTINE RETURNS TEMPERATURE (KELVIN) AND SPECIFIC HUMIDITY
C   (G/G) GIVEN VIRTUAL TEMPERATURE (KELVIN), PRESSURE (MB), AND
C   RELATIVE HUMIDITY (RANGE 0.0-1.0).

C   IT USES AN ITERATIVE NEWTON-RAPHSON TECHNIQUE.  FOUR ITERATIONS ARE
C   GENERALLY ADEQUATE TO PROVIDE CONVERGENCE TO 5 DECIMAL PLACES.
C   THE WOBUS FUNCTION FOR SATURATION VAPOR PRESSURE OVER LIQUID WATER
C   IS USED

C        S.G. BENJAMIN  13-DEC-1990        Original version
C        S.G. BENJAMIN  02-JAN-1995        New version

        IMPLICIT NONE
        REAL TV,RH,P,T,Q
        REAL TG,QG,FAC,FAC1,TVG,DTVDT,TGNEW
        REAL ESW,E
        INTEGER IT

        IT = 0
C -- TG - guess at non-virtual temperature
        E = RH*ESW(TV)
C   0.37803 = 0.6078 * 0.62197
        TG = TV/(1.+0.37803*E/(P-0.37803*E))
10      CONTINUE
        IT = IT + 1
C -- QG - guess at specific humidity
        E = RH*ESW(TG)
        QG = 0.62197*E/(P-0.37803*E)
        FAC = (1.+0.6078*QG)
C -- TVG - guess at virtual temperature
        TVG = TG*FAC
        E = RH*ESW(TG+0.1)
        FAC1= (1.+0.37803*E/(P-0.37803*E))
C -- DTVDT - d(Tv)/dT estimated over 0.1 K interval
C       from 1st term Taylor series expansion
        DTVDT = FAC + TG*(FAC1-FAC)/0.1
        TGNEW = TG + (TV-TVG)/DTVDT
C        PRINT *,',TV,TVG,TGNEW,QG,FAC,FAC1,DTVDT',TV,TVG,TGNEW,
C    1  QG,FAC,FAC1,DTVDT
        IF (ABS(TV-TVG).LT.0.01) GO TO 20
        IF (IT.GT.6) GO TO 20
        TG = TGNEW
        GO TO 10
20      CONTINUE
        T = TGNEW
        E = RH*ESW(T)
        Q = 0.62197*E/(P-0.37803*E)

        RETURN
        END



C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    FNDLVL      FIND THE CLOSEST ANALYSIS LEVEL TO OBS
C   PRGMMR: TRACY SMITH      ORG: FSL/PROFS  DATE: 90-07-19
C
C ABSTRACT: FIND WHAT ANALYSIS LEVEL IS CLOSEST TO THE OBSERVATION
C   FOR SINGLE LEVEL DATA.
C
C PROGRAM HISTORY LOG:
C
C
C USAGE:    CALL FNDLVL(OB,VCOORD,NZ,LVLC)
C   INPUT ARGUMENT LIST:
C     OB       - REAL    OBSERVATION RELATING TO VERTICAL COORDINATES
C     VCOORD   - REAL ARRAY  VERTICAL COORDINATES (LEVELS)
C     NZ       - INTEGER     NUMBER OF LEVELS
C
C   OUTPUT ARGUMENT LIST:
C     LVLC     - INTEGER     LEVEL THE OBSERVATION IS AT
C
C REMARKS:
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN-77
C   MACHINE:  NAS-9000
C$$$
      SUBROUTINE FNDLVL(OB,VERT_COORD,NZ,LVLC)

      IMPLICIT NONE

c     INCLUDE 'MAPSCON'

      REAL VERT_COORD(*)
C
      REAL OB,DIFF1,DIFF2
      INTEGER LVLC,I,NZ
C

C*** Decide if the parameter increases or decreases
        IF(VERT_COORD(1).LT.VERT_COORD(2)) THEN

C*** If the parameter decreases
      DO 70 I=1, NZ-1
          IF(OB.LT.VERT_COORD(I)) GO TO 71
 70   CONTINUE
 71   IF(I.NE.1) THEN
          DIFF1=VERT_COORD(I)-OB
          DIFF2=OB-VERT_COORD(I-1)
          IF(DIFF1.GT.DIFF2) THEN
             LVLC=I-1
          ELSE
             LVLC=I
          END IF
      ELSE
          LVLC=1
      END IF

C      WRITE(6,915) OB, LVLC, VERT_COORD(LVLC)
C 915  FORMAT(' CLOSEST LEVEL TO OBS=',F8.0,' IS ',I4,
C     ,        ' WHERE OBS=', F8.0)

C*** If the variable increases
        ELSE

      DO 80 I=1, NZ-1
          IF(OB.GT.VERT_COORD(I)) GO TO 81
 80   CONTINUE
 81   IF(I.NE.1) THEN
          DIFF1=OB-VERT_COORD(I)
          DIFF2=VERT_COORD(I-1)-OB
          IF(DIFF1.GT.DIFF2) THEN
             LVLC=I-1
          ELSE
             LVLC=I
          END IF
      ELSE
          LVLC=1
      END IF

C       WRITE(6,916) OB, LVLC, VERT_COORD(LVLC)
C916    FORMAT(' CLOSEST LEVEL TO OBS=',F8.0,' IS ',I4,
C     ,        ' WHERE OBS=', F8.0)
        ENDIF
C               ! the variable increases or decreases

      RETURN
      END



C$$$  SUBPROGRAM DOCUMENTATION BLOCK 
C                .      .    .                                       . 
C SUBPROGRAM:    SMOOTH      SMOOTH A METEOROLOGICAL FIELD
C   PRGMMR: STAN BENJAMIN    ORG: FSL/PROFS  DATE: 90-06-15 
C 
C ABSTRACT: SHAPIRO SMOOTHER. 
C 
C PROGRAM HISTORY LOG: 
C   85-12-09  S. BENJAMIN   ORIGINAL VERSION
C 
C USAGE:    CALL SMOOTH (FIELD,HOLD,IX,IY,SMTH) 
C   INPUT ARGUMENT LIST: 
C     FIELD    - REAL ARRAY  FIELD(IX,IY)
C                            METEOROLOGICAL FIELD
C     HOLD     - REAL ARRAY  HOLD(IX,2)
C                            HOLDING THE VALUE FOR FIELD
C     IX       - INTEGER     X COORDINATES OF FIELD
C     IY       - INTEGER     Y COORDINATES OF FIELD
C     SMTH     - REAL      
C
C   OUTPUT ARGUMENT LIST:   
C     FIELD    - REAL ARRAY  FIELD(IX,IY)
C                            SMOOTHED METEOROLOGICAL FIELD
C 
C REMARKS: REFERENCE: SHAPIRO, 1970: "SMOOTHING, FILTERING, AND
C   BOUNDARY EFFECTS", REV. GEOPHYS. SP. PHYS., 359-387.
C   THIS FILTER IS OF THE TYPE 
C         Z(I) = (1-S)Z(I) + S(Z(I+1)+Z(I-1))/2
C   FOR A FILTER WHICH IS SUPPOSED TO DAMP 2DX WAVES COMPLETELY
C   BUT LEAVE 4DX AND LONGER WITH LITTLE DAMPING,
C   IT SHOULD BE RUN WITH 2 PASSES USING SMTH (OR S) OF 0.5
C   AND -0.5.
C   
C ATTRIBUTES: 
C   LANGUAGE: FORTRAN-77 + EXTENSIONS
C   MACHINE:  NAS-9000, VAX, UNIX
C$$$ 
C**********************************************************************
C**********************************************************************

      SUBROUTINE SMOOTH (FIELD,HOLD,IX,IY,SMTH)

C**********************************************************************
C**********************************************************************

      REAL      FIELD(IX,IY),
     1  HOLD (IX,2)
      SMTH1 = 0.25 * SMTH * SMTH
      SMTH2 = 0.5  * SMTH * (1.-SMTH)
      SMTH3 = (1.-SMTH) * (1.-SMTH)
      SMTH4 = (1.-SMTH)
      SMTH5 = 0.5 * SMTH
      I1 = 2
      I2 = 1
       DO 30 J=2,IY-1
        IT = I1
        I1 = I2
        I2 = IT
         DO 10 I = 2,IX-1
           SUM1 = FIELD (I-1,J+1) + FIELD (I-1,J-1)
     1   + FIELD (I+1,J+1) + FIELD (I+1,J-1)
          SUM2 = FIELD (I  ,J+1) + FIELD (I+1,J  )
     1   + FIELD (I  ,J-1) + FIELD (I-1,J  )
          HOLD(I,I1) = SMTH1*SUM1 + SMTH2*SUM2 + SMTH3*FIELD(I,J)
10      CONTINUE
        IF (J.EQ.2) GO TO 200
         DO 20 I=2,IX-1
          FIELD(I,J-1) = HOLD(I,I2)
20      CONTINUE
200     CONTINUE
30     CONTINUE


       DO 40 I = 2,IX-1
        FIELD (I,IY-1) = HOLD(I,I1)
40    CONTINUE

       DO 50 I = 2,IX-1
        FIELD(I,1) = SMTH4* FIELD(I,1) 
     1       + SMTH5 * (FIELD(I-1,1) + FIELD(I+1,1))
        FIELD(I,IY) = SMTH4* FIELD(I,IY) 
     1       + SMTH5 * (FIELD(I-1,IY) + FIELD(I+1,IY))
50    CONTINUE

       DO 60 J = 2,IY-1
        FIELD(1,J) = SMTH4* FIELD(1,J) 
     1       + SMTH5 * (FIELD(1,J-1) + FIELD(1,J+1))
        FIELD(IX,J) = SMTH4* FIELD(IX,J) 
     1       + SMTH5 * (FIELD(IX,J-1) + FIELD(IX,J+1))
60      CONTINUE

      RETURN
       END


C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C SUBPROGRAM:    ESW         COMPUTE THE SATURATION VAPOR PRESSURE
C   PRGMMR: STAN BENJAMIN    ORG: FSL/PROFS  DATE: 90-07-19
C
C ABSTRACT: FUNCTION TO RETURNS THE SATURATION VAPOR PRESSURE
C   ESW (MILLIBARS) OVER LIQUID WATER GIVEN THE TEMPERATURE
C   T (CELSIUS OR KELVIN).
C
C PROGRAM HISTORY LOG:
C   82-05-17  TOM SCHLATTER   ORIGINAL VERSION
C
C USAGE:    ESW(T)
C   INPUT ARGUMENT LIST:
C     T        - REAL    TEMPERATURE (CELSIUS OR KELVIN)
C
C   OUTPUT ARGUMENT LIST:
C
C REMARKS: THE POLYNOMIAL APPROXIMATION BELOW IS DUE TO HERMAN WOBUS,
C   A MATHEMATICIAN WHO WORKED AT THE NAVY WEATHER RESEARCH FACILITY,
C   NORFOLK, VIRGINIA. THE COEFFICIENTS OF
C   THE POLYNOMIAL WERE CHOSEN TO FIT THE VALUES IN TABLE 94 ON
C   PP. 351-353 OF THE SMITHSONIAN METEOROLOGICAL TABLES BY ROLAND
C   LIST (6TH EDITION). THE APPROXIMATION IS VALID FOR
C   -50 < T < 100C.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN-77
C   MACHINE:  NAS-9000
C$$$
        FUNCTION ESW(T1)
C
C   THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE ESW (MILLIBARS)
C   OVER LIQUID WATER GIVEN THE TEMPERATURE T (CELSIUS).
C
C       BAKER,SCHLATTER 17-MAY-1982     Original version
C
C   THE POLYNOMIAL APPROXIMATION BELOW IS DUE TO HERMAN WOBUS, A 
C   MATHEMATICIAN WHO WORKED AT THE NAVY WEATHER RESEARCH FACILITY, 
C   NORFOLK, VIRGINIA. THE COEFFICIENTS OF THE 
C   POLYNOMIAL WERE CHOSEN TO FIT THE VALUES IN TABLE 94 ON PP. 351-353
C 
C   OF THE SMITHSONIAN METEOROLOGICAL TABLES BY ROLAND LIST (6TH 
C   EDITION). THE APPROXIMATION IS VALID FOR -50 < T < 100C.
C
C   ES0 = SATURATION VAPOR RESSURE OVER LIQUID WATER AT 0C
        DATA ES0/6.1078/
        T = T1
        IF(T.GT.150.) T=T-273.15
        POL = 0.99999683       + T*(-0.90826951E-02 +
     1     T*(0.78736169E-04   + T*(-0.61117958E-06 +
     2     T*(0.43884187E-08   + T*(-0.29883885E-10 +
     3     T*(0.21874425E-12   + T*(-0.17892321E-14 +
     4     T*(0.11112018E-16   + T*(-0.30994571E-19)))))))))
        ESW = ES0/POL**8
        RETURN
         END
