!copyright (C) 2001  MSC-RPN COMM  %%%RPNPHY%%%
***FONCTION STHTAW  -  CALCULE TW OU THETAW
*

      FUNCTION STHTAW(HU,TX,PX,LNPS,MODP,SWTT,SWPH,SWTH) 1,2
*
#include "impnone.cdk"
      REAL STHTAW, HU, TX, PX, LNPS
      INTEGER MODP
      LOGICAL SWTT, SWPH, SWTH
*
*Author
*          N. Brunet  (Jan91)
*Revision
* 001      B. Bilodeau  (August 1991)- Adaptation to UNIX
* 002      N. Brunet  (May 1994) - change to l2ocprv in order
*          to obtain results closer to those of tephigram
* 003      M. Lepine (March 2003) -  CVMG... Replacements
*
*Object
*          to return TW or THETAW by calculating from specific
*          humidity, temperature and pressure. Result returned is in
*          Kelvins
*
*Arguments
*
*          - Input -
* HU       specific humidity in kg/kg
* TX       temperature or virtual temperature in Kelvins
* PX       see MODP
* LNPS     see MODP
* MODP     pressure mode (SI units only):
*          =0; pressure level in PX
*          =1; sigma level in PX
*          =2; sigma level in PX and logarithm of sigma level in
*          LNPS
*          =3; all points of pressure in LNPS(NI,*) in Pascals
*          =4; sigma level in PX and logarithm of sigma level in
*          LNPS(in millibars unless using SI units)
*          =5; logarithm of pressure level in PX(in millibars unless
*          using SI units)
* SWTT     .TRUE. to pass TT for argument
*          .FALSE. to pass TV for argument
* SWPH     .TRUE. to consider water and ice phase
*          .FALSE. to consider water phase only
* SWTH     .TRUE. to calculate THETA
*          .FALSE. to calculate TW
*
*
*IMPLICITES
#include "consphy.cdk"
*MODULES
      EXTERNAL INCTPHY
*
**
*--------------------------------------------------------------------
      REAL Q1, DQ1, Q0, PN, TH, H, FT0, DFT0
      REAL D, DLP, QP, PR, FAC, L2OCPRV
      INTEGER ITER, N, J
*
#include "dintern.cdk"
#include "fintern.cdk"
*--------------------------------------------------------------------
#include "initcph.cdk"
*
      TH = TX
      IF(.NOT.SWTT)TH = FOTTV(TX,HU)
#include "modpr1.cdk"
*
*     TROUVE D'ABORD  TW
*
      Q0 = HU
*
*     SOLUTIONNE PAR METHODE DE NEWTON DU 1ER ORDRE.
*
      DO 10 ITER=1,4
         H = HTVOCP(TH)
            Q1=FOQST(TH,PN)
            IF(.NOT.SWPH)Q1=FOQSA(TH,PN)
            DQ1=FODQS(Q1,TH)
            IF(.NOT.SWPH)DQ1=FODQA(Q1,TH)
         FT0 = -H * (Q1-Q0)
         DFT0 = -1.0 - (H*DQ1)
         TH = TH - FT0/DFT0
         Q0 = Q0 + FT0/(DFT0*H)
10    CONTINUE
*
      IF(.NOT.SWTH)THEN
         STHTAW = TH
         RETURN
      END IF
*
*     TROUVE THETAW EN UTILISANT LA PENTE DE L'ADIABATIQUE
*     MOUILLEE ET EN FAISANT LES CALCULS PAR TRANCHE DE
*     1000 PA (PASSANT DE PN A 100000 PA).
*
      IF(PN.EQ.100000.)THEN
         STHTAW = TH
      ELSE
         DLP=1000.
         L2OCPRV = 1.35E+07
         D = 100000. - PN
         N = INT(ABS(D/DLP))
         IF(D.LT.0.)DLP = -DLP
         D = D - FLOAT(N)*DLP
         IF(D.NE.0.)THEN
            N = N + 1
         ELSE
            D = DLP
         END IF
*
         PR = PN
*
         DO 100 J=1,N
            QP = D/PR
            Q1 = FOQST(TH,PR)
            IF(.NOT.SWPH)Q1 = FOQSA(TH,PR)
            H = HTVOCP(TH)
            FAC = (CAPPA + (H*Q1/TH))/(1. + L2OCPRV*Q1/(TH*TH))
            TH = TH*(1. + QP*FAC)
*
            PR = PR + D
            IF(J.EQ.1)D = DLP
*
100      CONTINUE
*
         STHTAW = TH
*
      END IF
*
      CONTAINS
#include "flatscp.cdk"
      END