!copyright (C) 2002  MSC-RPN COMM  %%%RPNPHY%%%
***S/P LIN_ADJTQ
*
#include "phy_macros_f.h"

      SUBROUTINE LIN_ADJTQ  ( T, Q, P, NI ) 2
#include "impnone.cdk"
*
C
      INTEGER NI
      REAL T(NI), Q(NI), P(NI)
*
*Author
*          J.-F. Mahfouf (Sept 2002) 
*
*Revision
*
*Object
*          To calculate wet bulb temperature 
*          and associated specific humidity at saturation 
*          for an air parcel characterized by T and Q at pressure P
*          Newton method with two iterations
*
*          For an unsaturated parcel : final values = initial values
*
*Arguments
*
*            - Outputs/inputs -
* T        Temperature (K)
* Q        Specific Humidity  (kg/kg)
*
*            - Inputs -
* P        Pressure  (Pa)
* NI       Horizontal Dimension
*
C
      LOGICAL LO
      INTEGER jl
      REAL ZA1, CHLS, DELTA2
*
************************************************************************
*     AUTOMATIC ARRAYS
************************************************************************
*
      AUTOMATIC ( LO1    , LOGICAL , (NI   ))
*
      AUTOMATIC ( ZA3    , REAL    , (NI   ))
      AUTOMATIC ( ZA4    , REAL    , (NI   ))
      AUTOMATIC ( ZLDCP  , REAL    , (NI   ))
      AUTOMATIC ( ZT_0   , REAL    , (NI   ))
      AUTOMATIC ( ZT_1   , REAL    , (NI   ))
      AUTOMATIC ( ZQ_0   , REAL    , (NI   ))
      AUTOMATIC ( ZQ_1   , REAL    , (NI   ))
      AUTOMATIC ( ZESC_0 , REAL    , (NI   ))
      AUTOMATIC ( ZESC_1 , REAL    , (NI   ))    
      AUTOMATIC ( ZQSC_0 , REAL    , (NI   ))
      AUTOMATIC ( ZQSC_1 , REAL    , (NI   ))
      AUTOMATIC ( ZDESC_0, REAL    , (NI   ))
      AUTOMATIC ( ZDESC_1, REAL    , (NI   ))
      AUTOMATIC ( ZCOR_0 , REAL    , (NI   ))
      AUTOMATIC ( ZCOR_1 , REAL    , (NI   ))
      AUTOMATIC ( ZQCD_0 , REAL    , (NI   ))
      AUTOMATIC ( ZQCD_1 , REAL    , (NI   ))
*
************************************************************************

C
C*    PHYSICAL CONSTANTS.
C     -------- ----------
C
#include "consphy.cdk"
#include "dintern.cdk"
#include "fintern.cdk"
C
      DELTA2 = CPV/CPD - 1.
      CHLS   = CHLC + CHLF
      ZA1 = 610.78
C
      DO jl=1,NI
         IF ( T(jl) > TRPL ) THEN
            ZA3(jl)=17.269
            ZA4(jl)=35.860
            ZLDCP(jl) = CHLC / (CPD*(1.+DELTA2*Q(jl)))
         ELSE
            ZA3(jl)=21.875
            ZA4(jl)= 7.660
            ZLDCP(jl) = CHLS / (CPD*(1.+DELTA2*Q(jl)))
         ENDIF
      END DO
C
C*    Set-up initial values for Newton iterations
C
      DO jl=1,NI
         ZT_0(jl) = T(jl)
         ZQ_0(jl) = Q(jl)
      END DO
C
C*    First iteration
C
      DO jl=1,NI
         ZESC_0(jl)  = ZA1*EXP(ZA3(jl)*(ZT_0(jl)-TRPL)/
     *                 (ZT_0(jl)-ZA4(jl)))
         ZQSC_0(jl)  = EPS1*ZESC_0(jl)/(P(jl)-EPS2*ZESC_0(jl))
         ZDESC_0(jl) = ZA3(jl)*(TRPL-ZA4(jl))*ZESC_0(jl)/
     *                 ((ZT_0(jl)-ZA4(jl))**2)
         ZCOR_0(jl)  = ZLDCP(jl)*EPS1*P(jl)*ZDESC_0(jl)/
     *                 ((P(jl)-EPS2*ZESC_0(jl))**2)
         ZQCD_0(jl)  = (ZQ_0(jl)-ZQSC_0(jl))/(1.0+ZCOR_0(jl))
         IF ( ZQ_0(jl) < ZQSC_0(jl) ) THEN
            ZQCD_0(jl) = 0.0
         ENDIF
         ZQ_1(jl) = ZQ_0(jl) - ZQCD_0(jl)
         ZT_1(jl) = ZT_0(jl) + ZQCD_0(jl)*ZLDCP(jl)
         LO1(jl) = ZQCD_0(jl) /= 0.0
      END DO
C
      LO=.FALSE.
      DO jl=1,NI
         LO=LO.OR.LO1(jl)
      END DO
C
      IF (.NOT.LO) RETURN
C
C*    Second iteration
C
      DO jl=1,NI
         ZESC_1(jl)  = ZA1*EXP(ZA3(jl)*(ZT_1(jl)-TRPL)/
     *                 (ZT_1(jl)-ZA4(jl)))
         ZQSC_1(jl)  = EPS1*ZESC_1(jl)/(P(jl)-EPS2*ZESC_1(jl))
         ZDESC_1(jl) = ZA3(jl)*(TRPL-ZA4(jl))*ZESC_1(jl)/
     *                 ((ZT_1(jl)-ZA4(jl))**2)
         ZCOR_1(jl)  = ZLDCP(jl)*EPS1*P(jl)*ZDESC_1(jl)/
     *                 ((P(jl)-EPS2*ZESC_1(jl))**2)
         ZQCD_1(jl)  = (ZQ_1(jl)-ZQSC_1(jl))/(1.0+ZCOR_1(jl))
         IF ( ZQ_0(jl) < ZQSC_0(jl) ) THEN
            ZQCD_1(jl) = 0.0
         ENDIF
         Q(jl) = ZQ_1(jl) - ZQCD_1(jl)
         T(jl) = ZT_1(jl) + ZQCD_1(jl)*ZLDCP(jl)
      END DO
C      
      RETURN
      END SUBROUTINE LIN_ADJTQ