      PROGRAM INITGS9R                                                                                                              
C     PROGRAM INITGS9R(ICTL,       GPINIT,       GSINIT,       OUTPUT,  )       I2                                                  
C    1           TAPE1=ICTL, TAPE2=GPINIT, TAPE3=GSINIT, TAPE6=OUTPUT)                                                              
C     ----------------------------------------------------------------          I2                                                  
C                                                                               I2                                                  
C     MAY 04/12 - S.KHARIN (REVERSE VERT. INTERPOL. AND HYBRIDIZATION FOR Q)    I2                                                  
C     JUN 07/10 - S.KHARIN (MODIFY QMIN AS IN GSHUMH. MODIFY SLOPE NEAR SURF.)  I2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       I2                                                  
C     JUN 13/02 - M. LAZARE  (PREVIOUS ROUTINE INITGS8)                         I2                                                  
C                                                                               I2                                                  
CINITGS9R - INTERPOLATES SIGMA LEVEL GRIDS FROM PRESS. LEVELS FOR GCM14+ 2  1   I1                                                  
C                                                                               I3                                                  
CAUTHOR  - M.LAZARE                                                             I3                                                  
C                                                                               I3                                                  
CPURPOSE - CONVERTS A PRESSURE LEVEL GAUSSIAN GRID DATASET OF INITIAL           I3                                                  
C          MODEL VARIABLES (SFPR,PHI,TEMP,U,V,Q) TO A SIGMA/HYBRID LEVEL        I3                                                  
C          DATASET (LNSP,TEMP,U,V,ES) SUITABLE FOR SPECTRAL ANALYSIS.           I3                                                  
C          NOTE - IT CALLS PTSK WHICH CONVERTS SPEC. HUM. TO MODEL MOISTURE     I3                                                  
C                 VARIABLE 'MOIST' WHICH CONTROLS THE CHOICE OF MOISTURE        I3                                                  
C                 VARIABLE... 4HT-TD/4H  TD/4H   Q/4HSQRT/4HRLNQ  FOR           I3                                                  
C                 DEW POINT DEPRESSION/DEW POINT/SPEC.HUM./Q**0.5/-LN(Q)**-1    I3                                                  
C                                                                               I3                                                  
CINPUT FILES...                                                                 I3                                                  
C                                                                               I3                                                  
C      ICTL   = INITIALIZATION CONTROL DATASET (SEE ICNTRLH).                   I3                                                  
C      GPINIT = PRESSURE LEVEL GAUSSIAN GRIDS. ALL EXCEPT PS MUST               I3                                                  
C               BE AVAILABLE ON THE SAME PRESSURE LEVELS. NOTE THAT             I3                                                  
C               U,V MUST BE SCALED (MODEL) WIND COMPONENTS.                     I3                                                  
C               (I.E. (U,V)*COS(LAT)/A, WHERE A = EARTH RADIUS).                I3                                                  
C                                                                               I3                                                  
C                 NAME   VARIABLE          UNITS                                I3                                                  
C                                                                               I3                                                  
C                   PS   SURFACE PRESSURE  PA                                   I3                                                  
C                 TEMP   TEMPERATURE       DEG K                                I3                                                  
C                    U   BIG WIND COMPONENT    1./SEC                           I3                                                  
C                    V   BIG WIND COMPONENT    1./SEC                           I3                                                  
C                   ES   SPECIFIC HUMIDITY G/G                                  I3                                                  
C                                                                               I3                                                  
COUTPUT FILE...                                                                 I3                                                  
C                                                                               I3                                                  
C      GSINIT = SIGMA/HYBRID LEVEL GAUSSIAN GRIDS OF...                         I3                                                  
C                                                                               I3                                                  
C                 NAME   VARIABLE          UNITS     LEVELS                     I3                                                  
C                                                                               I3                                                  
C                 LNSP   LN(SF.PRES.)      PA           1   SFC                 I3                                                  
C                 TEMP   TEMPERATURE       DEG K     ILEV   SH                  I3                                                  
C                    U   WIND COMPONENT   1./SEC     ILEV   SG                  I3                                                  
C                    V   WIND COMPONENT   1./SEC     ILEV   SG                  I3                                                  
C                   ES   MOISTURE VARIABLE           ILEV   SH                  I3                                                  
C----------------------------------------------------------------------------                                                       
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
                                                                                                                                    
      PARAMETER (MAXL=500)                                                                                                          
                                                                                                                                    
      INTEGER LP(MAXL), LG(MAXL), LH(MAXL)                                                                                          
                                                                                                                                    
      LOGICAL OK                                                                                                                    
      CHARACTER*4 STRING                                                                                                            
                                                                                                                                    
      REAL, ALLOCATABLE ::                                                                                                          
     +     U(:),V(:),T(:),ES(:),GZS(:),SP(:)                                                                                        
                                                                                                                                    
      REAL, ALLOCATABLE ::                                                                                                          
     +     AG(:),BG(:),AH(:),BH(:),SG(:),                                                                                           
     +     SH(:),SIG(:),PR(:),PRLOG(:),                                                                                             
     +     FPCOL(:),DFDLNP(:),DLNP(:)                                                                                               
                                                                                                                                    
      INTEGER KBUF(16)                                                                                                              
      INTEGER, ALLOCATABLE :: IBUF(:)                                                                                               
      COMMON /MACHTYP/ MACHINE,INTSIZE                                                                                              
C-----------------------------------------------------------------------                                                            
      NFIL=4                                                                                                                        
      CALL JCLPNT(NFIL,1,2,3,6)                                                                                                     
                                                                                                                                    
C     * GET THE PRESSURE LEVEL VALUES AND THE ANALYSIS GRID SIZE                                                                    
C     * FROM THE TEMPERATURE FIELDS IN FILE GPINIT.                                                                                 
                                                                                                                                    
      REWIND 2                                                                                                                      
      CALL FIND(2,NC4TO8("GRID"),-1,NC4TO8("TEMP"),-1,OK)                                                                           
      IF(.NOT.OK) CALL                             XIT('INITGS9R',-1)                                                               
                                                                                                                                    
C     * GET ARRAY DIMENSIONS                                                                                                        
                                                                                                                                    
      CALL FBUFFIN(-2,KBUF,-8,K,LEN)                                                                                                
      BACKSPACE 2                                                                                                                   
      CALL PRTLAB (KBUF)                                                                                                            
      IJ=KBUF(5)*KBUF(6)*MACHINE                                                                                                    
C                                                                                                                                   
C     * DETERMINE NUMBER OF PRESSURE LEVELS FROM INPUT FILES                                                                        
C                                                                                                                                   
      CALL FILEV(LP,NPL,KBUF,-2)                                                                                                    
      IF((NPL.LT.1).OR.(NPL.GT.MAXL)) CALL         XIT('INITGS9R',-2)                                                               
      NPL1=NPL+1                                                                                                                    
C                                                                                                                                   
C     * DETERMINE NUMBER OF MODEL LEVELS FROM INPUT CONTROL FILE                                                                    
C                                                                                                                                   
      REWIND 1                                                                                                                      
      READ(1,END=903) LABL,NSL                                                                                                      
      IF((NSL.LT.1).OR.(NSL.GT.MAXL)) CALL         XIT('INITGS9R',-3)                                                               
      IJK=IJ*MAX(NPL,NSL)                                                                                                           
C                                                                                                                                   
C     * ALLOCATE ARRAYS                                                                                                             
C                                                                                                                                   
      ALLOCATE(IBUF(IJ+8*MACHINE),                                                                                                  
     +     AG(NSL),BG(NSL),AH(NSL),BH(NSL),SG(NSL),                                                                                 
     +     SH(NSL),SIG(NSL),PR(NPL),PRLOG(NPL),                                                                                     
     +     FPCOL(NPL),DFDLNP(NPL+1),DLNP(NPL),                                                                                      
     +     U(IJK),V(IJK),T(IJK),ES(IJK),GZS(IJ),SP(IJ))                                                                             
                                                                                                                                    
      CALL LVDCODE(PR,LP,NPL)                                                                                                       
      ILG =KBUF(5)-1                                                                                                                
      ILAT=KBUF(6)                                                                                                                  
      WRITE(6,6005) NPL,(PR(L),L=1,NPL)                                                                                             
                                                                                                                                    
C     * READ MODEL SIGMA VALUES FROM THE CONTROL FILE.                                                                              
                                                                                                                                    
      REWIND 1                                                                                                                      
      READ(1,END=903) LABL,NSL,(SG(L),L=1,NSL)                                                                                      
     1                        ,(SH(L),L=1,NSL),                                                                                     
     2                LAY,ICOORD,PTOIT,MOIST,SREF,SPOW                                                                              
      CALL WRITLEV(SG,NSL,' SG ')                                                                                                   
      CALL WRITLEV(SH,NSL,' SH ')                                                                                                   
      WRITE(6,6016) LAY,ICOORD,PTOIT                                                                                                
      WRITE(6,6017) MOIST,SREF,SPOW                                                                                                 
C                                                                                                                                   
      CALL SIGLAB2 (LG,LH, SG,SH, NSL)                                                                                              
      CALL COORDAB (AG,BG, NSL,SG, ICOORD,PTOIT)                                                                                    
      CALL COORDAB (AH,BH, NSL,SH, ICOORD,PTOIT)                                                                                    
                                                                                                                                    
C     * STOP IF THERE IS NOT ENOUGH SPACE.                                                                                          
                                                                                                                                    
      LEN = ILAT*(ILG+1)                                                                                                            
      IF(    LEN.GT.IJ ) CALL                      XIT('INITGS9R',-4)                                                               
      IF(NPL*LEN.GT.IJK) CALL                      XIT('INITGS9R',-5)                                                               
      IF(NSL*LEN.GT.IJK) CALL                      XIT('INITGS9R',-6)                                                               
                                                                                                                                    
C     * PERFORM THE VERTICAL INTERPOLATION, WITH CONVERSION OF ES.                                                                  
                                                                                                                                    
      CALL PTSJ2  (1,2,3,T,U,V,ES,GZS,SP,LEN,NPL,PR,NSL,                                                                            
     1             AG,BG,AH,BH,NPL+1,FPCOL,DFDLNP,DLNP,                                                                             
     2             LP,LG,LH,SIG,PRLOG,IBUF,IJ+8,MOIST,SREF,SPOW)                                                                    
                                                                                                                                    
      DEALLOCATE(IBUF,                                                                                                              
     +     AG,BG,AH,BH,SG,                                                                                                          
     +     SH,SIG,PR,PRLOG,                                                                                                         
     +     FPCOL,DFDLNP,DLNP,                                                                                                       
     +     U,V,T,ES,GZS,SP)                                                                                                         
      CALL                                         XIT('INITGS9R',0)                                                                
                                                                                                                                    
C     * E.O.F. ON FILE ICTL.                                                                                                        
                                                                                                                                    
  903 CALL                                         XIT('INITGS9R',-7)                                                               
C-----------------------------------------------------------------------                                                            
 6005 FORMAT('0  PRES LEVELS',I5,20F5.0/('0',18X,20F5.0))                                                                           
 6010 FORMAT(I6,A4,/,(10X,10F6.3))                                                                                                  
 6016 FORMAT(' LAY=',I5,', COORD=',A4,', P.LID=',F10.3,' (PA)')                                                                     
 6017 FORMAT(' MOIST=',A4,', SREF=',E12.4,', SPOW=',F10.3)                                                                          
      END                                                                                                                           
      SUBROUTINE PTSJ2                                                                                                              
     1           (NFC,NF1,NF2,T,U,V,ES,GZS,SP,LA,NPL,PR,NSL,                                                                        
     2            AG,BG,AH,BH,NPL1,FPCOL,DFDLNP ,DLNP ,                                                                             
     3            LP,LG,LH,SIG,PRLOG, IBUF,MAXP8,MOIST,SREF,SPOW)                                                                   
C                                                                                                                                   
C     * MAY 04/12 - S.KHARIN. - REVERSE VERT. INTERPOL. AND HYBRIDIZATION FOR Q.                                                    
C     * JUN 07/10 - S.KHARIN. - MODIFY QMIN FOR ES AS IN GSHUMH.                                                                    
C     * JUN 12/02 - M.LAZARE. - PASS SREF,SPOW AND USE IN QHYB,SLQB.                                                                
C     * OCT 19/01 - M.LAZARE. - PREVIOUS VERSION PTSI.                                                                              
C     * AUG 13/90 - M.LAZARE. - PREVIOUS VERSION PTSH.                                                                              
C                                                                                                                                   
C     * CONVERTS NPL PRESSURE LEVELS OF GAUSSIAN GRIDS (ILG1,ILAT)                                                                  
C     * ON FILE NF1 TO NSL ETA LEVELS ON FILE NF2 FOR HYBRID MODEL.                                                                 
C     * CONVERT ES=SHUM TO MODEL MOISTURE VARIABLE, CONTROLLED BY MOIST.                                                            
C     * ETAG FOR MOMENTUM MID LAYER, ETAH FOR THERMODYNAMIC MID LAYER.                                                              
C     * PRESSURE LEVELS IN PR ARE ORDERED TOP DOWN.                                                                                 
C     * (U,T),(V,ES) MAY BE EQUIVALENCED.                                                                                           
                                                                                                                                    
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
                                                                                                                                    
      REAL T(LA,NPL), U(LA,NPL), V(LA,NPL), ES(LA,NPL), GZS(LA), SP(LA)                                                             
                                                                                                                                    
      REAL AG    (NSL),BG    (NSL),AH     (NSL) ,BH   (NSL)                                                                         
      REAL PR    (NPL ),FPCOL (NPL ),DFDLNP (NPL1),DLNP (NPL)                                                                       
      REAL SIG   (NSL ),PRLOG (NPL )                                                                                                
      INTEGER LP (NPL ),LG    (NSL ),LH     (NSL )                                                                                  
      INTEGER IBUF(MAXP8)                                                                                                           
                                                                                                                                    
      COMMON/EPS / A, B, EPS1, EPS2                                                                                                 
      COMMON/EPSICE/ AICE,BICE,TICE,QMIN                                                                                            
      COMMON/HTCP  / T1S,T2S,AI,BI,AW,BW,SLP                                                                                        
C                                                                                                                                   
C     * STATEMENT FUNCTIONS.                                                                                                        
C                                                                                                                                   
C     * A) COMPUTES THE FRACTIONAL PROBABILITY OF WATER PHASE                                                                       
C     *    EXISTING AS A FUNCTION OF TEMPERATURE (FROM ROCKEL,                                                                      
C     *    RASCHKE AND WEYRES, BEITR. PHYS. ATMOSPH., 1991.)                                                                        
                                                                                                                                    
      FRACW(TTT) = MERGE( 1.E0,                                                                                                     
     &  0.0059E0+0.9941E0*EXP(-0.003102E0*(T1S-TTT)**2),                                                                            
     &  TTT.GE.T1S )                                                                                                                
                                                                                                                                    
C     * B) COMPUTES THE DEW-POINT TEMPERATURE OVER WATER OR ICE.                                                                    
                                                                                                                                    
      TDW(XXX) = B/(A - XXX)                                                                                                        
      TDI(XXX) = BICE/(AICE - XXX)                                                                                                  
      TDEFF(XXX,TTT) = FRACW(TTT)*TDW(XXX) + (1.E0-FRACW(TTT))*TDI(XXX)                                                             
C--------------------------------------------------------------------                                                               
      MAXX = MAXP8-8                                                                                                                
                                                                                                                                    
C     * GET THERMODYNAMIC CONSTANTS.                                                                                                
                                                                                                                                    
      CALL SPWCON7(FVORT,PI)                                                                                                        
                                                                                                                                    
C     * SET NPL TO PRESSURE VALUES. SET PRLOG TO LN(PR).                                                                            
                                                                                                                                    
      DO 110 L=1,NPL                                                                                                                
         PRLOG(L) = LOG(PR(L))                                                                                                      
  110 CONTINUE                                                                                                                      
C---------------------------------------------------------------------                                                              
C     * READ MOUNTAIN GEOPOTENTIAL (M/SEC)**2.                                                                                      
C                                                                                                                                   
C      READ(NFC)                                                                                                                    
C      INSTEAD USE RDSKIP TO SKIP RECORDS.                                                                                          
      CALL RDSKIP(NFC,NSKIP,OK)                                                                                                     
      CALL GETFLD2(NFC,GZS,NC4TO8("GRID"),-1,NC4TO8("PHIS"),1,                                                                      
     1                                           IBUF,MAXX,OK)                                                                      
      IF(.NOT.OK) CALL                             XIT('PTSJ2',-1)                                                                  
      CALL PRTLAB(IBUF)                                                                                                             
      REWIND NF1                                                                                                                    
      NREC=0                                                                                                                        
C                                                                                                                                   
C     * READ SURFACE PRESSURE AND CONVERT TO LSNP.                                                                                  
C                                                                                                                                   
  200 CALL GETFLD2(NF1,SP,NC4TO8("GRID"),-1,NC4TO8("  PS"),1,                                                                       
     1                                          IBUF,MAXX,OK)                                                                       
      IF(.NOT.OK) THEN                                                                                                              
        IF(NREC.EQ.0)CALL                          XIT('PTSJ2',-2)                                                                  
        RETURN                                                                                                                      
      ENDIF                                                                                                                         
      CALL PRTLAB(IBUF)                                                                                                             
      DO 210 I=1,LA                                                                                                                 
         SP(I)=LOG(SP(I))                                                                                                           
  210 CONTINUE                                                                                                                      
                                                                                                                                    
C     * READ TEMPERATURES (DEG K).                                                                                                  
                                                                                                                                    
      DO 220 L=1,NPL                                                                                                                
         CALL GETFLD2(NF1, T(1,L),NC4TO8("GRID"),-1,NC4TO8("TEMP"),                                                                 
     1                                         LP (L),IBUF,MAXX,OK)                                                                 
         IF(.NOT.OK) CALL                          XIT('PTSJ2',-3)                                                                  
         CALL PRTLAB(IBUF)                                                                                                          
  220 CONTINUE                                                                                                                      
                                                                                                                                    
C     * FROM DT/DP=(R*T)/(P*CP),  R/CP=2./7.,  AND T=280.                                                                           
C     * WE GET RLDN = DT/D(LN P) = 80.                                                                                              
C     * LTES = (0,1) TO PUT TEMP,ES ON (HALF,FULL) LEVELS.                                                                          
                                                                                                                                    
      RLUP =  0.E0                                                                                                                  
      RLDN = 80.E0                                                                                                                  
      CALL PAEL2 (T,LA,SIG,NSL, T,PRLOG,NPL,SP,RLUP,RLDN,                                                                           
     1            AH,BH,NPL+1,FPCOL,DFDLNP ,DLNP )                                                                                  
                                                                                                                                    
C     * OUTPUT PHIS,LNSP AND TEMPERATURES.                                                                                          
                                                                                                                                    
      IBUF(4)=1                                                                                                                     
      IBUF(3)=NC4TO8("PHIS")                                                                                                        
      CALL PUTFLD2(NF2, GZS,IBUF,MAXX)                                                                                              
      CALL PRTLAB(IBUF)                                                                                                             
      IBUF(3)=NC4TO8("LNSP")                                                                                                        
      CALL PUTFLD2(NF2, SP ,IBUF,MAXX)                                                                                              
      CALL PRTLAB(IBUF)                                                                                                             
      DO 280 K=1,NSL                                                                                                                
         IBUF(3)=NC4TO8("TEMP")                                                                                                     
         IBUF(4)=LH(K)                                                                                                              
         CALL PUTFLD2(NF2,  T(1,K),IBUF,MAXX)                                                                                       
         CALL PRTLAB(IBUF)                                                                                                          
  280 CONTINUE                                                                                                                      
                                                                                                                                    
C---------------------------------------------------------------------                                                              
C     * READ MODEL WINDS (U,V)*COS(LAT)/(EARTH RADIUS).                                                                             
                                                                                                                                    
      DO 510 L=1,NPL                                                                                                                
         CALL GETFLD2(NF1, U(1,L),NC4TO8("GRID"),-1,NC4TO8("   U"),                                                                 
     1                                         LP (L),IBUF,MAXX,OK)                                                                 
         IF(.NOT.OK) CALL                          XIT('PTSJ2',-4)                                                                  
         CALL PRTLAB(IBUF)                                                                                                          
  510 CONTINUE                                                                                                                      
      DO 520 L=1,NPL                                                                                                                
         CALL GETFLD2(NF1, V(1,L),NC4TO8("GRID"),-1,NC4TO8("   V"),                                                                 
     1                                         LP (L),IBUF,MAXX,OK)                                                                 
         IF(.NOT.OK) CALL                          XIT('PTSJ2',-5)                                                                  
         CALL PRTLAB(IBUF)                                                                                                          
  520 CONTINUE                                                                                                                      
                                                                                                                                    
C     * INTERPOLATE TO ETA LEVELS.                                                                                                  
                                                                                                                                    
      RLUP = 0.E0                                                                                                                   
      RLDN = 0.E0                                                                                                                   
      CALL PAEL2( U,LA, SIG,NSL, U,PRLOG,NPL,SP,RLUP,RLDN,                                                                          
     1            AG,BG,NPL+1,FPCOL,DFDLNP ,DLNP )                                                                                  
      CALL PAEL2( V,LA, SIG,NSL, V,PRLOG,NPL,SP,RLUP,RLDN,                                                                          
     1            AG,BG,NPL+1,FPCOL,DFDLNP ,DLNP )                                                                                  
                                                                                                                                    
C     * OUTPUT PAIRS OF U,V ON FILE NF2.                                                                                            
                                                                                                                                    
      DO 580 K=1,NSL                                                                                                                
         IBUF(4)=LG(K)                                                                                                              
                                                                                                                                    
         IBUF(3)=NC4TO8("   U")                                                                                                     
         CALL PUTFLD2(NF2,  U(1,K),IBUF,MAXX)                                                                                       
         CALL PRTLAB(IBUF)                                                                                                          
                                                                                                                                    
         IBUF(3)=NC4TO8("   V")                                                                                                     
         CALL PUTFLD2(NF2,  V(1,K),IBUF,MAXX)                                                                                       
         CALL PRTLAB(IBUF)                                                                                                          
  580 CONTINUE                                                                                                                      
                                                                                                                                    
C---------------------------------------------------------------------                                                              
C     * READ SPECIFIC HUMIDITY (G/G) FROM FILE NF1.                                                                                 
C     * LTES = (0,1) TO PUT TEMP,ES ON (HALF,FULL) LEVELS.                                                                          
                                                                                                                                    
      DO 410 L=1,NPL                                                                                                                
         CALL GETFLD2(NF1,ES(1,L),NC4TO8("GRID"),-1,NC4TO8("SHUM"),                                                                 
     1                                         LP (L),IBUF,MAXX,OK)                                                                 
         IF(.NOT.OK) CALL                          XIT('PTSJ2',-6)                                                                  
         CALL PRTLAB(IBUF)                                                                                                          
  410 CONTINUE                                                                                                                      
                                                                                                                                    
C     * CONVERT TO MODEL MOISTURE VARIABLE.                                                                                         
                                                                                                                                    
      DO 470 N=1,NPL                                                                                                                
       IF(MOIST.EQ.NC4TO8("RLNQ")) THEN                                                                                             
          DO 430 I=1,LA                                                                                                             
             SHUM   = ES(I,N)                                                                                                       
             ES(I,N)= -1.E0/LOG(SHUM)                                                                                               
  430     CONTINUE                                                                                                                  
                                                                                                                                    
       ELSEIF(MOIST.EQ.NC4TO8("SQRT")) THEN                                                                                         
          DO 435 I=1,LA                                                                                                             
             SHUM   = ES(I,N)                                                                                                       
             ES(I,N)= SQRT(SHUM)                                                                                                    
  435     CONTINUE                                                                                                                  
                                                                                                                                    
       ELSEIF(MOIST.EQ.NC4TO8("   Q") .OR. MOIST.EQ.NC4TO8("SL3D")) THEN                                                            
C                                                                                                                                   
C         * NO CONVERSION REQUIRED.                                                                                                 
C                                                                                                                                   
                                                                                                                                    
       ELSEIF(MOIST.EQ.NC4TO8("  TD")) THEN                                                                                         
          DO 445 I=1,LA                                                                                                             
             PSPA   = 100.E0*EXP(SP(I))                                                                                             
             CALL NIVCAL (SIGMA,AH(N),BH(N),PSPA,1,1,1)                                                                             
             PRES   = 0.01E0*PSPA*SIGMA                                                                                             
             ALVP   = LOG(ES(I,N)*PRES/(EPS1+EPS2*ES(I,N)))                                                                         
             TD     = TDEFF(ALVP,T(I,N))                                                                                            
             ES(I,N)= TD                                                                                                            
  445     CONTINUE                                                                                                                  
                                                                                                                                    
       ELSEIF(MOIST.EQ.NC4TO8(" LNQ")) THEN                                                                                         
          DO 450 I=1,LA                                                                                                             
             ES(I,N)= LOG(SHUM)                                                                                                     
  450     CONTINUE                                                                                                                  
                                                                                                                                    
       ELSEIF(MOIST.EQ.NC4TO8("QHYB") .OR. MOIST.EQ.NC4TO8("SLQB")) THEN                                                            
          PINV=1.E0/SPOW                                                                                                            
C SK: ADJUST QMIN TO MAKE IT CONSISTENT WITH THAT IN GSHUMH                                                                         
          QMIN=MAX(SREF*1.E-16,1.E-36)                                                                                              
          DO 460 I=1,LA                                                                                                             
             SHUM   = MAX(ES(I,N),QMIN)                                                                                             
             IF(SHUM.GE.SREF) THEN                                                                                                  
                ES(I,N)= SHUM                                                                                                       
             ELSE                                                                                                                   
                ES(I,N)= SREF/((1.E0+SPOW*LOG(SREF/SHUM))**PINV)                                                                    
             ENDIF                                                                                                                  
  460     CONTINUE                                                                                                                  
                                                                                                                                    
       ELSEIF(MOIST.EQ.NC4TO8("T-TD")) THEN                                                                                         
          DO 465 I=1,LA                                                                                                             
             PSPA   = 100.E0*EXP(SP(I))                                                                                             
             CALL NIVCAL (SIGMA,AH(N),BH(N),PSPA,1,1,1)                                                                             
             PRES   = 0.01E0*PSPA*SIGMA                                                                                             
             ALVP   = LOG(ES(I,N)*PRES/(EPS1+EPS2*ES(I,N)))                                                                         
             TD     = TDEFF(ALVP,T(I,N))                                                                                            
             ES(I,N)= T(I,N)-TD                                                                                                     
  465     CONTINUE                                                                                                                  
                                                                                                                                    
       ELSE                                                                                                                         
         CALL                                      XIT('PTSJ2',-50)                                                                 
       ENDIF                                                                                                                        
                                                                                                                                    
  470 CONTINUE                                                                                                                      
                                                                                                                                    
C     * INTERPOLATE TO ETA LEVELS. ES MUST BE POSITIVE EVERYWHERE                                                                   
C     * (BOUND IS 3 PPMV).                                                                                                          
                                                                                                                                    
      RLUP = 0.E0                                                                                                                   
      RLDN = 0.E0                                                                                                                   
      CALL PAEL2 (ES,LA,SIG,NSL, ES,PRLOG,NPL,SP,RLUP,RLDN,                                                                         
     1            AH,BH,NPL+1,FPCOL,DFDLNP ,DLNP )                                                                                  
      DO 425 N=1,NSL                                                                                                                
      DO 425 I=1,LA                                                                                                                 
C SK: ALREADY HYBRIDIZED                                                                                                            
CC       ES(I,N)=MAX(ES(I,N),2.E-6)                                                                                                 
         PSPA   = 100.E0*EXP(SP(I))                                                                                                 
         CALL NIVCAL (SIGMA,AH(N),BH(N),PSPA,1,1,1)                                                                                 
         PRES   = 0.01E0*PSPA*SIGMA                                                                                                 
  425 CONTINUE                                                                                                                      
                                                                                                                                    
C     * OUTPUT ES ON FILE NF2.                                                                                                      
                                                                                                                                    
      DO 480 K=1,NSL                                                                                                                
         IBUF(3)=NC4TO8("  ES")                                                                                                     
         IBUF(4)=LH(K)                                                                                                              
         CALL PUTFLD2(NF2, ES(1,K),IBUF,MAXX)                                                                                       
         CALL PRTLAB(IBUF)                                                                                                          
  480 CONTINUE                                                                                                                      
                                                                                                                                    
      NREC=NREC+1                                                                                                                   
      GOTO 200                                                                                                                      
C--------------------------------------------------------------------                                                               
 6026 FORMAT(' ',60X,A4,I6,1X,A4,5I10)                                                                                              
      END                                                                                                                           
      SUBROUTINE PAEL2(FS, LA,SIG,NSL, FP,PRLOG,NPL, PSLOG,RLUP,RLDN,                                                               
     1                 A,B,NPL1,FPC,DFDLNP,DLNP)                                                                                    
                                                                                                                                    
C     * DERIVE RLDN FROM THE BOTTOM TWO LAYERS - S.KHARIN                                                                           
C     * FEB 02/88 - R.LAPRISE.                                                                                                      
C     * INTERPOLATES MULTI-LEVEL SET OF GRIDS FROM PRESSURE LEVELS                                                                  
C     * TO ETA LEVELS FOR HYBRID MODEL.                                                                                             
C     * THE INTERPOLATION IS LINEAR IN LN(PRES).                                                                                    
C                                                                                                                                   
C     * ALL GRIDS HAVE THE SAME HORIZONTAL SIZE (LA POINTS).                                                                        
C     * FP    = INPUT  GRIDS ON PRESSURE LEVELS.                                                                                    
C     * FS    = OUTPUT GRIDS ON ETA      LEVELS.                                                                                    
C     *     (NOTE THAT FP AND FS MAY BE EQUIVALENCED IN CALLING PROGRAM)                                                            
C     * PSLOG = INPUT  GRID OF LN(SURFACE PRESSURE IN MB).                                                                          
C                                                                                                                                   
C     * PRLOG(NPL) = VALUES OF INPUT LOG(PRESSURE LEVEL (MB)).                                                                      
C     * SIG(NSL) = WORKING ARRAY FOR PRESSURE ON ETA LEVELS OF OUTPUT FIELD.                                                        
C     *            (BOTH MUST BE MONOTONIC AND INCREASING).                                                                         
C     * NPL1  = NPL+1.                                                                                                              
C                                                                                                                                   
C     * RLUP  = LAPSE RATE USED TO EXTRAPOLATE ABOVE PRLOG(1).                                                                      
C     * RLDN  = LAPSE RATE USED TO EXTRAPOLATE BELOW PRLOG(NSL).                                                                    
C     *         UNITS OF RLUP,RLDN ARE  DF/D(LN PRES).                                                                              
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      REAL FP   (LA,NPL),FS (LA,NSL),PSLOG (LA)                                                                                     
      REAL PRLOG(   NPL)                                                                                                            
                                                                                                                                    
      REAL A(NSL),B(NSL)                                                                                                            
                                                                                                                                    
C     * WORK SPACE.                                                                                                                 
                                                                                                                                    
      REAL FPC(NPL),DFDLNP(NPL1),DLNP(NPL),SIG(NSL)                                                                                 
C---------------------------------------------------------------------                                                              
C     * PRECOMPUTE INVERSE OF DELTA LN PRES.                                                                                        
C                                                                                                                                   
      DO 180 L=1,NPL-1                                                                                                              
         DLNP(L) =1.E0 / (PRLOG(L+1)-PRLOG(L))                                                                                      
  180 CONTINUE                                                                                                                      
                                                                                                                                    
C     * LOOP OVER ALL HORIZONTAL POINTS.                                                                                            
                                                                                                                                    
      DO 500 I=1,LA                                                                                                                 
                                                                                                                                    
C     * GET A COLUMN OF FP (ON PRESSURE LEVELS).                                                                                    
                                                                                                                                    
      DO 210 L=1,NPL                                                                                                                
         FPC(L)=FP(I,L)                                                                                                             
  210 CONTINUE                                                                                                                      
                                                                                                                                    
C     * COMPUTE VERTICAL DERIVATIVE OVER ALL PRESSURE INTERVALS.                                                                    
                                                                                                                                    
      DO 260 L=1,NPL-1                                                                                                              
  260 DFDLNP(L+1  ) = (FPC(L+1)-FPC(L)) *DLNP(L)                                                                                    
      DFDLNP(1    ) = RLUP                                                                                                          
      DFDLNP(NPL+1) = RLDN                                                                                                          
                                                                                                                                    
C     * COMPUTE LOG(PRESSURE(MB)) ON MODEL LEVELS                                                                                   
                                                                                                                                    
      DO 300 L=1,NSL                                                                                                                
         SIG(L)=LOG(A(L)/100.E0+B(L)*EXP(PSLOG(I)))                                                                                 
  300 CONTINUE                                                                                                                      
                                                                                                                                    
C     * LOOP OVER SIGMA LEVELS TO BE INTERPOLATED.                                                                                  
                                                                                                                                    
      K=1                                                                                                                           
      DO 350 N=1,NSL                                                                                                                
                                                                                                                                    
C     * FIND WHICH PRESSURE INTERVAL WE ARE IN.                                                                                     
                                                                                                                                    
      DO 310 L=K,NPL                                                                                                                
      INTVL=L                                                                                                                       
  310 IF(SIG(N).LT.PRLOG(L)) GO TO 320                                                                                              
      INTVL=NPL+1                                                                                                                   
  320 K=INTVL-1                                                                                                                     
      IF(K.EQ.0) K=1                                                                                                                
                                                                                                                                    
C     * NOW INTERPOLATE AT THIS POINT.                                                                                              
C     * SK: USE THE GRADIENT IN THE ABOVE LAYER TO EXTRAPOLATE BELOW THE LOWEST LEVEL                                               
                                                                                                                                    
      FS(I,N)  = FPC(K)+DFDLNP(MIN(INTVL,NPL))*(SIG(N)-PRLOG(K))                                                                    
  350 CONTINUE                                                                                                                      
                                                                                                                                    
  500 CONTINUE                                                                                                                      
                                                                                                                                    
      RETURN                                                                                                                        
      END                                                                                                                           
