!copyright (C) 2001  MSC-RPN COMM  %%%RPNPHY%%%
*
*   DEFINITION DES FONCTIONS THERMODYNAMIQUES DE BASE
*   POUR LES CONSTANTES, UTILISER LE COMMON /CTESDYN/
*     NOTE: TOUTES LES FONCTIONS TRAVAILLENT AVEC LES UNITES S.I.
*           I.E. TTT EN DEG K, PRS EN PA, QQQ EN KG/KG
*          *** N. BRUNET - MAI 90 ***
*          * REVISION 01 - MAI 94 - N. BRUNET
*                          NOUVELLE VERSION POUR FAIBLES PRESSIONS
*          * REVISION 02 - AOUT 2000 - J-P TOVIESSI
*                          CALCUL EN REAL*8
*          * REVISION 03 - SEPT 2000 - N. BRUNET
*                          AJOUT DE NOUVELLES FONCTIONS
*          * REVISION 04 - JANV 2000 - J. MAILHOT
*                          FONCTIONS EN PHASE MIXTE
*          * REVISION 05 - DEC 2001 - G. LEMAY
*                          DOUBLE PRECISION POUR PHASE MIXTE
*          * REVISION 06 - AVR 2002 - A. PLANTE
*                          AJOUT DES NOUVELLES FONCTIONS FOTTVH ET FOTVHT
*
*     FONCTION DE TENSION DE VAPEUR SATURANTE (TETENS) - EW OU EI SELON TT
      FOEW(TTT) = 610.78D0*DEXP( DMIN1(DSIGN(17.269D0,
     W DBLE(TTT)-DBLE(TRPL)),DSIGN
     W (21.875D0,DBLE(TTT)-DBLE(TRPL)))*DABS(DBLE(TTT)-DBLE(TRPL))/
     W (DBLE(TTT)-35.86D0+DMAX1(0.D0,DSIGN
     W (28.2D0,DBLE(TRPL)-DBLE(TTT)))))
*
*     FONCTION CALCULANT LA DERIVEE SELON T DE  LN EW (OU LN EI)
      FODLE(TTT)=(4097.93D0+DMAX1(0.D0,DSIGN(1709.88D0,
     W DBLE(TRPL)-DBLE(TTT))))
     W /((DBLE(TTT)-35.86D0+DMAX1(0.D0,DSIGN(28.2D0,
     W DBLE(TRPL)-DBLE(TTT))))*(DBLE(TTT)-35.86D0+DMAX1(0.D0
     W ,DSIGN(28.2D0,DBLE(TRPL)-DBLE(TTT)))))
*
*     FONCTION CALCULANT L'HUMIDITE SPECIFIQUE SATURANTE (QSAT)
      FOQST(TTT,PRS) = DBLE(EPS1)/(DMAX1(1.D0,DBLE(PRS)/FOEW(TTT))-
     W DBLE(EPS2))
*
*     FONCTION CALCULANT LA DERIVEE DE QSAT SELON T
      FODQS(QST,TTT)=DBLE(QST)*(1.D0+DBLE(DELTA)*DBLE(QST))*FODLE(TTT)
*     QST EST LA SORTIE DE FOQST
*
*     FONCTION CALCULANT TENSION VAP (EEE) FN DE HUM SP (QQQ) ET PRS
      FOEFQ(QQQ,PRS) = DMIN1(DBLE(PRS),(DBLE(QQQ)*DBLE(PRS)) /
     W (DBLE(EPS1) + DBLE(EPS2)*DBLE(QQQ)))
*
*      FONCTION CALCULANT HUM SP (QQQ) DE TENS. VAP (EEE) ET PRES (PRS)
      FOQFE(EEE,PRS) = DMIN1(1.D0,DBLE(EPS1)*DBLE(EEE)/(DBLE(PRS)-
     W DBLE(EPS2)*DBLE(EEE)))
*
*      FONCTION CALCULANT TEMP VIRT. (TVI) DE TEMP (TTT) ET HUM SP (QQQ)
      FOTVT(TTT,QQQ) = DBLE(TTT) * (1.0D0 + DBLE(DELTA)*DBLE(QQQ))

*      FONCTION CALCULANT TEMP VIRT. (TVI) DE TEMP (TTT), HUM SP (QQQ) ET
*      MASSE SP DES HYDROMETEORES.
      FOTVHT(TTT,QQQ,QQH) = DBLE(TTT) * 
     +     (1.0D0 + DBLE(DELTA)*DBLE(QQQ) - DBLE(QQH))
*
*      FONCTION CALCULANT TTT DE TEMP VIRT. (TVI) ET HUM SP (QQQ)
      FOTTV(TVI,QQQ) = DBLE(TVI) / (1.0D0 + DBLE(DELTA)*DBLE(QQQ))

*      FONCTION CALCULANT TTT DE TEMP VIRT. (TVI), HUM SP (QQQ) ET
*      MASSE SP DES HYDROMETEORES (QQH)
      FOTTVH(TVI,QQQ,QQH) = DBLE(TVI) / 
     +     (1.0D0 + DBLE(DELTA)*DBLE(QQQ) - DBLE(QQH))
*
*      FONCTION CALCULANT HUM REL DE HUM SP (QQQ), TEMP (TTT) ET PRES (PRS)
*      HR = E/ESAT
       FOHR(QQQ,TTT,PRS) = DMIN1(DBLE(PRS),FOEFQ(QQQ,PRS)) / FOEW(TTT)
*
*     FONCTION CALCULANT LA CHALEUR LATENTE DE CONDENSATION
      FOLV(TTT) =DBLE(CHLC) - 2317.D0*(DBLE(TTT)-DBLE(TRPL))
*
*     FONCTION CALCULANT LA CHALEUR LATENTE DE SUBLIMATION
      FOLS(TTT) = DBLE(CHLC)+DBLE(CHLF)+(DBLE(CPV)-
     +            (7.24D0*DBLE(TTT)+128.4D0))*(DBLE(TTT)-DBLE(TRPL))
*
*     FONCTION RESOLVANT L'EQN. DE POISSON POUR LA TEMPERATURE
*     NOTE: SI PF=1000*100, "FOPOIT" DONNE LE THETA STANDARD
      FOPOIT(T00,PR0,PF)=DBLE(T00)*(DBLE(PR0)/DBLE(PF))**
     +                 (-DBLE(CAPPA))
*
*     FONCTION RESOLVANT L'EQN. DE POISSON POUR LA PRESSION
      FOPOIP(T00,TF,PR0)=DBLE(PR0)*DEXP(-(DLOG(DBLE(T00)/DBLE(TF))/
     +                 DBLE(CAPPA)))
*
*     LES 5 FONCTIONS SUIVANTES SONT VALIDES DANS LE CONTEXTE OU ON
*     NE DESIRE PAS TENIR COMPTE DE LA PHASE GLACE DANS LES CALCULS
*     DE SATURATION.
*   FONCTION DE VAPEUR SATURANTE (TETENS)
      FOEWA(TTT)=610.78D0*DEXP(17.269D0*(DBLE(TTT)-DBLE(TRPL))/
     W (DBLE(TTT)-35.86D0))
*   FONCTION CALCULANT LA DERIVEE SELON T DE LN EW
      FODLA(TTT)=17.269D0*(DBLE(TRPL)-35.86D0)/(DBLE(TTT)-35.86D0)**2
*   FONCTION CALCULANT L'HUMIDITE SPECIFIQUE SATURANTE
      FOQSA(TTT,PRS)=DBLE(EPS1)/(DMAX1(1.D0,DBLE(PRS)/FOEWA(TTT))-
     W DBLE(EPS2))
*   FONCTION CALCULANT LA DERIVEE DE QSAT SELON T
      FODQA(QST,TTT)=DBLE(QST)*(1.D0+DBLE(DELTA)*DBLE(QST))*FODLA(TTT)
*   FONCTION CALCULANT L'HUMIDITE RELATIVE
      FOHRA(QQQ,TTT,PRS)=DMIN1(DBLE(PRS),FOEFQ(QQQ,PRS))/FOEWA(TTT)
*
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
*
*   Definition of basic thermodynamic functions in mixed-phase mode
*     FFF is the fraction of ice and DDFF its derivative w/r to T
*     NOTE: S.I. units are used
*           i.e. TTT in deg K, PRS in Pa
*          *** J. Mailhot - Jan. 2000 ***
*
*     Saturation calculations in presence of liquid phase only
*     Function for saturation vapor pressure (TETENS)
      FESI(TTT)=610.78D0*DEXP(21.875D0*(DBLE(TTT)-DBLE(TRPL))/
     %       (DBLE(TTT)-7.66D0)  )
      FDLESI(TTT)=21.875D0*(DBLE(TRPL)-7.66D0)/(DBLE(TTT)-7.66D0)**2
      FESMX(TTT,FFF) = (1.D0-DBLE(FFF))*FOEWA(TTT)+DBLE(FFF)*FESI(TTT)
      FDLESMX(TTT,FFF,DDFF) = ( (1.D0-DBLE(FFF))*FOEWA(TTT)*FODLA(TTT)
     1                      + DBLE(FFF)*FESI(TTT)*FDLESI(TTT)
     1            + DBLE(DDFF)*(FESI(TTT)-FOEWA(TTT)) )/FESMX(TTT,FFF)
      FQSMX(TTT,PRS,FFF) = DBLE(EPS1)/
     %        (DMAX1(1.D0,DBLE(PRS)/FESMX(TTT,FFF) ) - DBLE(EPS2)  )
      FDQSMX(QSM,DLEMX) = DBLE(QSM ) *(1.D0 + DBLE(DELTA)* DBLE(QSM ) )
     %                     * DBLE(DLEMX )
*