      SUBROUTINE RELHUM4 (RH,T,Q,PCP,PRESSG,                            
     1                   PBLT,DEL,SHJ,                                  
     2                   ILEV,MILEV,ILG,IL1,IL2,MSG,                    
     3                   P,PCPDH,SPHJ,DQR,DQMAX,                        
     4                   HTHRESH,DTH,DQH,QSJ,DQJ,                       
     5                   HT,H,HS,EPH,Y,                                 
     6                   HZ,EJ,DQEX,DQEP         )                      

C     * AUG 16/95 - F. MAJAESS. REPLACE AMIN1/AMAX1 CALLS BY MIN/MAX.
C     * DEC 14/94 - M.LAZARE.  USE CONSTANT "L/CP" TO ENSURE BALANCING
C     *                        OF (BALT-BEG) OVER THE COURSE OF A YEAR.
C     * FEB 02/93 - M.LAZARE.  PREVIOUS CHANGE TAKEN OUT.
C     * DEC 18/92 - M.LAZARE.  HTHRESH NOW CALCULATED USING WATER/ICE
C     *                        CURVE.
C     * JUL 27/92 - M.LAZARE.  LIKE PREVIOUS VERSION RELHUM3 EXCEPT:                   
C     *                        A. USE OF "FRACW" ONSET CURVE TO DEFINE            
C     *                           WATER/ICE BOUNDARY, AS IN NEW
C     *                           CLOUDS5.
C     *                        B. RELATIVE HUMIDITY FIELD ("RH") DEFINED
C     *                           BEFORE RELEASE OF LARGE-SCALE
C     *                           CONDENSATION, AS IN GCM6.
                                                                        
C     * MAR 19/91 - M.LAZARE.  PREVIOUS VERSION RELHUM3 FOR GCM7U.                
                                                                        
C     * REMOVES SUPERSATURATION BASED ON LARGE SCALE EVAPORATION                  
C     * /CONDENSATION.                                                            
                                                                        
C     **************** INDEX OF VARIABLES ***********************                 
C     *  I    => INPUT ARRAYS.                                                    
C     * I/O   => INPUT/OUTPUT ARRAYS.                                             
C     *  W    => WORK ARRAYS.                                                     
C     * IC    => INPUT DATA CONSTANTS.                                            
C     *  C    => DATA CONSTANTS PERTAINING TO SUBROUTINE ITSELF.                  
                                                                        
C  I  * DEL      LOCAL SIGMA HALF-LEVEL THICKNESS (I.E. DSHJ).                    
C  W  * DQEX     MOISTURE REMOVED BY CONDENSATION.                                
C  W  * DQH      MOISTURE CHANGE CAUSED BY A RELEASE OF LATENT HEAT.              
C  W  * DQJ      SATURATION DEFICIT.                                              
C  W  * DQMAX    AMOUNT OF SUB-SATURATION USED TO DETERMINE EVAPORATION.          
C  W  * DTH      TEMPERATURE CHANGE CAUSED BY A RELEASE OF LATENT HEAT.           
C  W  * DQR      MOISTURE EQUIVALENT TO PRECIPITATION AMOUNT.                     
C  W  * EJ       VAPOUR PRESSURE.                                                 
C  W  * EPH      FACTOR TO CALCULATE MOISTURE RELEASE UPON CONDENSATION.          
C  W  * H        RELATIVE HUMIDITY.                                               
C  W  * HS       CRITICAL SATURATION RELATIVE HUMIDITY.                           
C  W  * HT       RATIO OF LATENT HEAT OF VAPOUR TO SPECIFIC HEAT.                 
C  W  * HTHRESH  THRESHOLD RH FOR LARGE-SCALE CONDENSATION.                       
C  W  * HZ       FACTOR TO CALCULATE MOISTURE RELEASE UPON CONDENSATION.          
C IC  * ILEV     NUMBER OF MODEL LAYERS.                                          
C IC  * ILG      LON+2 = SIZE OF GRID SLICE.                                      
C IC  * IL1,IL2  START AND END LONGITUDE INDICES FOR LATITUDE CIRCLE.             
C IC  * MSG      NUMBER OF MISSING MOISTURE LEVELS AT THE TOP OF MODEL.           
C  W  * P        GRID SLICE OF AMBIENT PRESSURE IN MBS.                           
C I   * PBLT     ROW OF TOP-OF-PBL INDICES.                                       
C I/O * PCP      ROW OF PRECIPITATION RATE.                                       
C  W  * PCPDH    SCALED SURFACE PRESSURE.                                         
C  I  * PRESSG   ROW OF SURFACE PRESSURE IN PA.                                   
C I/O * Q        GRID SLICE OF SPECIFIC HUMIDITY.                                 
C  W  * QSJ      CRITICAL SPECIFIC HUMIDITY.                                      
C  O  * RH       CALCULATED GRID SLICE OF TEMPERATURE.                            
C  I  * SHJ      GRID SLICE OF LOCAL HALF-LEVEL SIGMA VALUES.                     
C  W  * SPHJ     SPECIFIC HUMIDITY.                                               
C I/O * T        GRID SLICE OF TEMPERATURE.                                       
C  W  * Y        FACTOR TO CALCULATE MOISTURE RELEASE UPON CONDENSATION.          
                                                                        
C     ***********************************************************                 
                                                                        
C     * MULTI-LEVEL I/O FIELDS:                                                   
                                                                        
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      REAL RH(ILG,MILEV),   T(ILG,MILEV),    Q(ILG,MILEV),              
     1    DEL(ILG,MILEV), SHJ(ILG,MILEV)                                
                                                                        
C     * SINGLE-LEVEL I/O FIELDS:                                                  
                                                                        
      REAL     PCP(ILG),  PRESSG(ILG),    PBLT(ILG)                     
                                                                        
C     * WORK FIELDS:                                                              
                                                                        
      REAL P(ILG,MILEV)                                                 
                                                                        
      REAL   PCPDH(ILG),    SPHJ(ILG),     DQR(ILG),   DQMAX(ILG)       
      REAL HTHRESH(ILG),     DTH(ILG),     DQH(ILG),     QSJ(ILG)       
      REAL     DQJ(ILG),      HT(ILG),       H(ILG),      HS(ILG)       
      REAL     EPH(ILG),       Y(ILG),      HZ(ILG),      EJ(ILG)       
      REAL    DQEX(ILG),    DQEP(ILG)                                   
                                                                        
      COMMON/ADJPCP/ HC,HF,HM,AA,DEPTH,LHEAT,MOIADJ,MOIFLX              
      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 DEFENITIONS.                                          
                                                                        
C     * A) CRITICAL SATURATION RELATIVE HUMIDITY.                                 
                                                                        
      CR(HM,HHH,AA)     =     MERGE( AA*(HHH-HM)**3,                    
     &  AA*(2.E0-HM-HHH)**3+HHH-1.E0,        
     &  HHH.LE.1.E0)                         
      CRIRLH(HM,HHH,AA) = HHH -
     &  MERGE( CR(HM,HHH,AA), HHH-1.E0, HM+HHH.LE.2)
                                                                        
C     * B) 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     * C) COMPUTES THE RATIO OF LATENT HEAT OF VAPORIZATION OF                   
C     *    WATER OR ICE TO THE SPECIFIC HEAT OF AIR AT CONSTANT                   
C     *    PRESSURE CP. NOTE THAT A CONSTANT "L/CP" IS REQUIRED
C     *    TO MATCH SIMILAR VALUES USED ELSEWHERE IN THE CODE,
C     *    OTHERWISE A FICTITIOUS ENERGY SOURCE IS GENERATED WHICH
C     *    CAUSES (BALT-BEG) TO NOT BALANCE TO ZERO OVER THE COURSE
C     *    OF A YEAR.                                                           

      HTVOCP(TTT) = 2.501E6/1004.64E0                                   
C     TW(TTT)     = AW-BW*TTT                                                     
C     TI(TTT)     = AI-BI*TTT                                                     
C     HTVOCP(TTT) = FRACW(TTT)*TW(TTT) + (1.-FRACW(TTT))*TI(TTT)                  
                                                                        
C     * D) COMPUTES THE SATURATION VAPOUR PRESSURE OVER WATER OR ICE.             
                                                                        
      ESW(TTT)    = EXP(A-B/TTT)                                        
      ESI(TTT)    = EXP(AICE-BICE/TTT)                                  
      ESTEFF(TTT) = FRACW(TTT)*ESW(TTT) + (1.E0-FRACW(TTT))*ESI(TTT)    
                                                                        
C     * E) COMPUTES THE EFFECTIVE THRESHOLD FOR LARGE-SCALE 
C     *    CONDENSATION.
                                                                        
      HTHEFF(TTT,UUU,VVV) = FRACW(TTT)*UUU + (1.E0-FRACW(TTT))*VVV
C----------------------------------------------------------------------           
C     * INITIALIZE NECESSARY ARRAYS.                                              
C                                                                                 
      DO 50 IL=IL1,IL2                                                  
         PCP(IL)   = 0.E0                                               
         DQR(IL)   = 0.E0                                               
         PCPDH(IL) = PRESSG(IL)*DEPTH                                   
   50 CONTINUE                                                          
C                                                                                 
      DO 60 L=1,ILEV                                                    
      DO 60 IL=IL1,IL2                                                  
         P(IL,L)   = SHJ(IL,L) * PRESSG(IL) * 0.01E0                    
   60 CONTINUE                                                          
C                                                                                 
C     ****************************************************************            
C     * CONVECTIVE OR STABLE HEATING/COOLING BY LARGE-SCALE                       
C     * CONDENSATION/EVAPORATION.                                                 
C     ****************************************************************            
      DO 200 L=1,ILEV                                                   
          IF (L.GT.MSG)                                      THEN       
              DO 100 IL=IL1,IL2                                         
                  EJ(IL)      = ESTEFF(T(IL,L))                         
                  SPHJ(IL)    = EPS1*EJ(IL)/(P(IL,L)-EPS2*EJ(IL))       
                  H(IL)       = MAX(Q(IL,L),QMIN)/SPHJ(IL)              
                  RH(IL,L)    = MIN(H(IL),1.E0)                         
                  HTHRESH(IL) = HM
                  IF(H(IL).GT.HTHRESH(IL))                          THEN
                     HS(IL)   = CRIRLH(HTHRESH(IL),H(IL),AA)            
                  ELSE                                                  
                     HS(IL)   = MIN(H(IL),1.E0)                         
                  ENDIF                                                 
                  QSJ(IL)     = HS(IL) * SPHJ(IL)                       
                  DQJ(IL)     = QSJ(IL) - Q(IL,L)                       
                  IF(DQJ(IL).GT.0.E0) THEN
                     DQMAX(IL) = MAX(HTHRESH(IL)*SPHJ(IL)-Q(IL,L),0.E0) 
                     DQEP (IL) = MAX(-DQR(IL)/DEL(IL,L), 0.E0)          
                     DQJ  (IL) = MIN(DQMAX(IL),DQEP(IL))                
                  ENDIF                                                 
                  HT(IL)     = HTVOCP(T(IL,L))                          
                  EPH(IL)    = QSJ(IL)*EPS2/(HS(IL)*EPS1)               
                  Y(IL)      = B/(T(IL,L)*T(IL,L))*QSJ(IL)*
     &              (1.E0+EPH(IL)) 
                  HZ(IL)     = HT(IL)*Y(IL)/
     &              (T(IL,L)*(1.E0+HT(IL)*Y(IL))) 
                  DQH(IL)    = DQJ(IL)/(1.E0+HT(IL)*Y(IL))              
     1                        * (1.E0+DQJ(IL)/Y(IL)*HZ(IL)*HZ(IL)       
     2                        * (B*(EPH(IL)+0.5E0)-T(IL,L)))            
                  DTH(IL)    = -HT(IL)*DQH(IL)                          
                  T(IL,L)    = T(IL,L)+DTH(IL)                          
                  Q(IL,L)    = Q(IL,L)+DQH(IL)                          
C                                                                                 
C                 * REMOVES LEFT OVER SUPERSATURATION TO INSURE                   
C                 * MAX(R.H.)<=1.                                                 
C                                                                                 
                  HT(IL)     = HTVOCP(T(IL,L))                          
                  EJ(IL)     = ESTEFF(T(IL,L))                          
                  SPHJ(IL)   = EPS1*EJ(IL)/(P(IL,L)-EPS2*EJ(IL))        
                  DQEX(IL)   = MAX(Q(IL,L) - SPHJ(IL), 0.E0)            
                  DQH(IL)    = DQH(IL) - DQEX(IL)                       
                  Q(IL,L)    = Q(IL,L) - DQEX(IL)                       
                  T(IL,L)    = T(IL,L) + HT(IL) * DQEX(IL)              
C                                                                                 
C                 * CONVECTIVE OR STABLE PRECIPITATION.                           
C                                                                                 
                  DQR(IL) = DQR(IL) + DEL(IL,L) * DQH(IL)               
  100         CONTINUE                                                  
          ENDIF                                                         
  200 CONTINUE                                                          
C                                                                                 
C     * UPDATE PRECIPITATION.                                                     
C                                                                                 
      DO 300 IL=IL1,IL2                                                 
          PCP(IL) = PCP(IL) - PCPDH(IL)*DQR(IL)                         
  300 CONTINUE                                                          
                                                                        
      RETURN                                                            
      END
