      SUBROUTINE EATHLT (TEMPS,PHIS,RGASMS,LNSPS,LA,NSL,NSL1,TH,NTHL,   
     1                   RLUP,RLDN,A,B,TEMPTH,PHITH,PRESTH,DPDTH,
     2                   TEMPC,PHIC,PRESS,THS,RMEAN,GAMMA,NU)
                                                                        
C     * OCT  5/94 - J.KOSHYK (INTRODUCE COLUMN ARRAYS TEMPC/PHIC TO CORRECT    
C                   CALCULATIONS OF TEMPTH/PHITH).        
C     * JUL 27/94 - J. KOSHYK  
                                                  
C     * INTERPOLATES MULTI-LEVEL SETS OF TEMP, PHI, PRES, DPDTHETA ON ETA        
C     * OR PRESSURE LEVELS TO THETA LEVELS.                                      
C                                                                                 
C     * ALL GRIDS HAVE THE SAME HORIZONTAL SIZE (LA POINTS).                      
C     * TEMPS    = INPUT TEMPERATURE GRIDS ON ETA/PRES LEVELS.                  
C     * PHIS     = INPUT GEOPOTENTIAL HEIGHT GRIDS ON ETA/PRES LEVELS.            
C     * RGASMS   = INPUT MOIST GAS CONSTANT GRIDS ON ETA/PRES LEVELS.             
C     * LNSPS    = INPUT GRID OF LN(SURFACE PRESSURE IN MB).                    
C     * NSL      = NUMBER OF ETA/PRES LEVELS.                                 
C     * TH(NTHL) = VALUES OF INPUT THETA LEVELS (K);
C     *            (MUST BE MONOTONIC AND DECREASING).
C     * RLUP     = -DT/DZ USED TO EXTRAPOLATE ABOVE TOP ETA/PRES.                
C     * RLDN     = -DT/DZ USED TO EXTRAPOLATE BELOW BOTTOM ETA/PRES.
C     * A, B     = PARAMETERS OF ETA VERTICAL DISCRETIZATION.
C     * TEMPTH   = OUTPUT GRIDS OF TEMPERATURE ON THETA LEVELS.
C     * PHITH    = OUTPUT GRIDS OF GEOPOTENTIAL HEIGHT ON THETA LEVELS.
C     * PRESTH   = OUTPUT GRIDS OF PRESSURE ON THETA LEVELS.
C     * DPDTH    = OUTPUT GRIDS OF DPRES/DTHETA ON THETA LEVELS.
C                                                                                 
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      REAL TEMPS  (LA,NSL),PHIS  (LA,NSL),RGASMS (LA,NSL),LNSPS     (LA)
      REAL TH       (NTHL),A        (NSL),B         (NSL) 
      REAL TEMPTH(LA,NTHL),PHITH(LA,NTHL),PRESTH(LA,NTHL),DPDTH(LA,NTHL)
                                                                        
C     * WORK SPACE.                                                               
                                                                        
      REAL TEMPC(NSL),PHIC(NSL)
      REAL PRESS (NSL),THS   (NSL)
      REAL RMEAN(NSL1),GAMMA(NSL1),NU(NSL1),KAPPA

      PARAMETER(P0=101320.E0, KAPPA=0.285714E0, G=9.81E0)

C---------------------------------------------------------------------------------
C     * LOOP OVER ALL HORIZONTAL POINTS.                                          

      DO 500 I=1,LA                                                     
                                                                        
C     * COMPUTE PRESSURE AND POTENTIAL TEMPERATURE ON ETA LEVELS.

      DO 150 L=1,NSL
        TEMPC(L)  = TEMPS(I,L)
        PHIC(L)   = PHIS(I,L)
        PRESS(L)   = A(L)+B(L)*(100.E0*EXP(LNSPS(I)))
        THS(L)     = TEMPC(L)*(P0/PRESS(L))**KAPPA
  150 CONTINUE

C     * COMPUTE THE LAPSE RATE, GAMMA = -DT/D PHI OVER ALL INPUT INTERVALS,
C     * AND CALCULATE RELATED QUANTITIES FOR LATER USE.

      DO 180 L=1,NSL-1                                                  
        RMEAN(L+1) =  (RGASMS(I,L+1)+RGASMS(I,L))/2.E0
        GAMMA(L+1) = -(TEMPC(L+1)-TEMPC(L))/(PHIC(L+1)-PHIC(L))
        NU   (L+1) =  1.E0/(KAPPA-RMEAN(L+1)*GAMMA(L+1))
180   CONTINUE                                                          
                                                                        
C     * ASSIGN VALUES OF QUANTITIES ABOVE HIGHEST AND BELOW LOWEST
C     * ETA/PRES LEVEL.

      GAMMA(1)     = RLUP/G
      GAMMA(NSL1)  = RLDN/G

C     * ASSUME VALUE OF RGASM ABOVE HIGHEST LEVEL = VALUE WITHIN HIGHEST LAYER,
C     *    AND VALUE OF RGASM BELOW LOWEST  LEVEL = VALUE WITHIN LOWEST  LAYER.

      RMEAN(1)     = RMEAN(2)
      RMEAN(NSL1)  = RMEAN(NSL)
      NU   (1)     = 1.E0/(KAPPA-RMEAN(1)*GAMMA(1))
      NU   (NSL1)  = 1.E0/(KAPPA-RMEAN(NSL1)*GAMMA(NSL1))
                                                                        
C     * LOOP OVER THETA LEVELS TO BE INTERPOLATED.                             
                                                                        
      K=1                                                               
      DO 350 N=1,NTHL                                                   
                                                                        
C     * FIND WHICH SIGMA INTERVAL WE ARE IN.                              
                                                                        
      DO 310 L=K,NSL                                                    
      INTVL=L                                                           
  310 IF(TH(N).GT.THS(L)) GO TO 320                                     
      INTVL=NSL+1                                                       
  320 K=INTVL-1                                                         
      IF(K.EQ.0) K=1                                                    
                                                                        
C     * NOW INTERPOLATE AT THIS POINT.                                            
      
      PRESTH(I,N) = PRESS(K)*(THS(K)/TH(N))**NU(INTVL)
      DPDTH (I,N) = -NU(INTVL)*(PRESTH(I,N)/TH(N))
      TEMPTH(I,N) = TH(N)*(PRESTH(I,N)/P0)**KAPPA
     
      IF (GAMMA(INTVL) .NE. 0) THEN
        PHITH (I,N) = PHIC(K)-(TEMPTH(I,N)-TEMPC(K))/GAMMA(INTVL)
      ELSE

C       * COMPUTE PHI FROM HYDROSTATIC EQUATION, ASSUMING ISOTHERMAL LAYER.
 
        PHITH (I,N) = PHIC(K)-RMEAN(INTVL)*TEMPTH(I,N)
     1                        *LOG(PRESTH(I,N)/PRESS(K))
      ENDIF
  
  350 CONTINUE
  500 CONTINUE                                                          
                                                                        
      RETURN                                                            
      END  
