      SUBROUTINE TQTEND2(DQDT,DSDT,Q,S,QS,QU,SU,MC,DU,QHAT,SHAT,DP,      
     1                   MU,MD,SD,QD,QL,DQLDT,DSR,
     2                   DSUBCLD,MB,JT,MX,
     3                   DT,ILEV,ILG,IL1G,IL2G,CP,MSG             )                
C
C     * MAY 03/2018 - N.MCFARLANE. ADDED CODE TO PREVENT STATCLD CRASHES.
C     * MAY 27/2006 - M.LAZARE.    NEW VERSION FOR GCM15F:
C     *                            - NOW USES INTERNAL WORK ARRAYS
C     *                              INSTEAD OF PASSED WORKSPACE.    
C     * MAY 05/2003 - K.VONSALZEN. PREVIOUS VERSION TQTEND:
C     *                            - CORRECT PROBLEMS WITH MOISTURE AND
C     *                              ENERGY CONSERVATION (I.E. AIR DENSITY).
C     * SEP 21/2001 - K.VONSALZEN. PREVIOUS VERSION Q1Q26.
C
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
                                                                             
      REAL       DQDT(ILG,ILEV),   DSDT(ILG,ILEV)                              
      REAL      DQLDT(ILG,ILEV), DSR(ILG,ILEV)
      REAL          Q(ILG,ILEV),      S(ILG,ILEV),     QS(ILG,ILEV),            
     1             QU(ILG,ILEV),     SU(ILG,ILEV),     MC(ILG,ILEV),            
     2             DU(ILG,ILEV),   QHAT(ILG,ILEV),   SHAT(ILG,ILEV),
     3             DP(ILG,ILEV),     MU(ILG,ILEV),     MD(ILG,ILEV),
     4             SD(ILG,ILEV),     QD(ILG,ILEV),     QL(ILG,ILEV)  
                                                                            
      REAL    DSUBCLD(ILG),  MB(ILG)                  
                                                                              
      INTEGER  JT(ILG),    MX(ILG)
                                                                              
C     * INTERNAL WORK FIELDS:                                                          
                 
      REAL  , DIMENSION(ILG,ILEV)  :: FACQ,FACS,DSDS,DQDQ
      REAL  , DIMENSION(ILG)       :: RL, ACOR
C----------------------------------------------------------------------------  
      RNU=1.
      DO 20 K=MSG+1,ILEV                                                     
      DO 20 IL=IL1G,IL2G                                                     
        DSDT(IL,K)=0.                                                       
        DQDT(IL,K)=0.                                                        
        DQLDT(IL,K)=0. 
        IF(K.GE.JT(IL) .AND. K.LT.ILEV)                             THEN    
          SDIF=0.                                                            
          IF(S(IL,K) .GT. 0. .OR. S(IL,K+1) .GT. 0.)                         
     1      SDIF=ABS((S(IL,K)-S(IL,K+1))/MAX(S(IL,K),S(IL,K+1)))             
          IF(SDIF. GT. 1.E-3)                               THEN              
            DSDS(IL,K)=   S(IL,K+1)/(S(IL,K)-S(IL,K+1))                       
     1                 -LOG( S(IL,K)/S(IL,K+1) )                             
     2                 *( S(IL,K+1)/(S(IL,K)-S(IL,K+1)) )**2                  
          ELSE                                                                
            DSDS(IL,K)=0.5                                                    
          ENDIF                                                               
          QDIF=0.                                                             
          IF(Q(IL,K) .GT. 0. .OR. Q(IL,K+1) .GT. 0.)                          
     1      QDIF=ABS((Q(IL,K)-Q(IL,K+1))/MAX(Q(IL,K),Q(IL,K+1)))             
          IF(QDIF .GT. 1.E-3)                               THEN               
            DQDQ(IL,K)=   Q(IL,K+1)/(Q(IL,K)-Q(IL,K+1))                        
     1                 -LOG( Q(IL,K)/Q(IL,K+1) )                              
     2                 *( Q(IL,K+1)/(Q(IL,K)-Q(IL,K+1)) )**2                   
          ELSE                                                                 
            DQDQ(IL,K)=0.5                                                     
          ENDIF                                                                
          FACS(IL,K)=1./(  1.+MC(IL,K+1)*DT*DSR(IL,K)
     1                    /DP(IL,K)*RNU*DSDS(IL,K) )         
          FACQ(IL,K)=1./(  1.+MC(IL,K+1)*DT*DSR(IL,K)
     1                    /DP(IL,K)*RNU*DQDQ(IL,K) )         
        ENDIF                                                                  
   20 CONTINUE
C
      DO 100 J=MSG+1,ILEV-1
      DO 100 IL=IL1G,IL2G                                                     
        IF(J.EQ.JT(IL) .AND. MB(IL).GT.0.)                          THEN      
          RL (IL)=2.501E6                                                     
          DSDT(IL,J)=FACS(IL,J)*DSR(IL,J)/DP(IL,J)                           
     1         * (MU(IL,J+1)*(SU(IL,J+1)-SHAT(IL,J+1)              
     2                             -RL(IL)/CP*QL(IL,J+1)       ) +           
     3            MD(IL,J+1)*(SD(IL,J+1)-SHAT(IL,J+1)) )           
          DQDT(IL,J)=FACQ(IL,J)*DSR(IL,J)/DP(IL,J)                           
     1         * (MU(IL,J+1)*(QU(IL,J+1)-QHAT(IL,J+1)              
     2                             +QL(IL,J+1)                 ) +           
     3            MD(IL,J+1)*(QD(IL,J+1)-QHAT(IL,J+1)) )           
        ENDIF                                                                
  100 CONTINUE                                                                
C                   
      BETA=1.                                                                 
      DO 200 J=MSG+1,ILEV                                                      
      DO 200 IL=IL1G,IL2G                                                    
        IF(J.GT.JT(IL) .AND. J.LT.MX(IL) .AND. MB(IL).GT.0.)        THEN     
          RL (IL)=2.501E6                                                    
          SHAT(IL,J)=SHAT(IL,J)+DSDT(IL,J-1)*DT*RNU*DSDS(IL,J-1)             
          QHAT(IL,J)=QHAT(IL,J)+DQDT(IL,J-1)*DT*RNU*DQDQ(IL,J-1)             
          DSDT(IL,J)=FACS(IL,J)                                              
     1               *DSR(IL,J)*((MC(IL,J)*(SHAT(IL,J)-S(IL,J))                        
     2               +MC(IL,J+1)*(S(IL,J)-SHAT(IL,J+1)))/DP(IL,J)            
     3     -RL(IL)/CP*DU(IL,J)*(BETA*QL(IL,J)+(1-BETA)*QL(IL,J+1))  )        
          DQDT(IL,J)=FACQ(IL,J)                                              
     1               *DSR(IL,J)*((MC(IL,J)*(QHAT(IL,J)-Q(IL,J))                        
     2               +MC(IL,J+1)*(Q(IL,J)-QHAT(IL,J+1)))/DP(IL,J)            
     3               +DU(IL,J)*(QS(IL,J)-Q(IL,J))                            
     4               +DU(IL,J)*(BETA*QL(IL,J)+(1-BETA)*QL(IL,J+1)) )         
        ENDIF                                                                
  200 CONTINUE      
C
      DO 300 J=MSG+1,ILEV                                                    
      DO 300 IL=IL1G,IL2G                                                    
        IF(J.EQ.MX(IL) .AND. MB(IL).GT.0.)                          THEN     
          SHAT(IL,J)=SHAT(IL,J)+DSDT(IL,J-1)*DT*RNU*DSDS(IL,J-1)             
          QHAT(IL,J)=QHAT(IL,J)+DQDT(IL,J-1)*DT*RNU*DQDQ(IL,J-1)             
          DSDT(IL,J)=(1./DSUBCLD(IL))                                          
     1              *( MU(IL,J)*( SHAT(IL,J)-SU(IL,J))          
     2              +MD(IL,J)*( SHAT(IL,J)-SD(IL,J)) )          
          DQDT(IL,J)=(1./DSUBCLD(IL))                                          
     1              *( MU(IL,J)*( QHAT(IL,J)-QU(IL,J))          
     2              +MD(IL,J)*( QHAT(IL,J)-QD(IL,J)) )          
        ENDIF                                                                
  300 CONTINUE      
C
C     * Rescale tendencies to protect against pathological cases where small
C     * humidity values may be made negative because of undershooting negative
C     * humidity tendencies. The same scaling must be applied to both humidity
C     * and temperature tendencies to preserve conservation constraints
C
      DO 320 IL=IL1G,IL2G
         ACOR (IL)=1.
  320 CONTINUE
C
      DO 325 L=MSG+1,ILEV
      DO 325 IL=IL1G,IL2G
         QTEMP=Q(IL,L)+DT*DQDT(IL,L)
         IF ( QTEMP.LT.0. ) THEN
            ACORT=-Q(IL,L)/(DQDT(IL,L)*DT)
            ACOR1=MAX(MIN(ACOR(IL),ACORT),0.)
            ACOR(IL)=(3.-2.*ACOR1)*ACOR1**3
         ENDIF
  325 CONTINUE
C
      DO 330 L=MSG+1,ILEV
      DO 330 IL=IL1G,IL2G
         DSDT(IL,L)=DSDT(IL,L)*ACOR(IL)
         DQDT(IL,L)=DQDT(IL,L)*ACOR(IL)
  330 CONTINUE
C                                                                             
      RETURN                                                                  
      END
