c-------
c comsat 
c-------
c
c     Statement functions and associated parameters.
c
c----------------------------------------------------------------------
c
c        Polynomials for svp(t), d(svp)/dt over water and ice are from
c        Lowe(1977),JAM,16,101-103.
c
      PARAMETER (ASAT0 =  6.107799961,
     *           ASAT1 =  4.436518521E-1,
     *           ASAT2 =  1.428945805E-2,
     *           ASAT3 =  2.650648471E-4,
     *           ASAT4 =  3.031240396E-6,
     *           ASAT5 =  2.034080948E-8,
     *           ASAT6 =  6.136820929E-11 )
c
      PARAMETER (BSAT0 =  6.109177956,
     *           BSAT1 =  5.03469897E-1,
     *           BSAT2 =  1.886013408E-2,
     *           BSAT3 =  4.176223716E-4,
     *           BSAT4 =  5.824720280E-6,
     *           BSAT5 =  4.838803174E-8,
     *           BSAT6 =  1.838826904E-10 )
c
      PARAMETER (CSAT0 =  4.438099984E-1,
     *           CSAT1 =  2.857002636E-2,
     *           CSAT2 =  7.938054040E-4,
     *           CSAT3 =  1.215215065E-5,
     *           CSAT4 =  1.036561403E-7,
     *           CSAT5 =  3.532421810E-10,
     *           CSAT6 = -7.090244804E-13 )
c
      PARAMETER (DSAT0 =  5.030305237E-1,
     *           DSAT1 =  3.773255020E-2,
     *           DSAT2 =  1.267995369E-3,
     *           DSAT3 =  2.477563108E-5,
     *           DSAT4 =  3.005693132E-7,
     *           DSAT5 =  2.158542548E-9,
     *           DSAT6 =  7.131097725E-12 )
c
c----------
c
c         Statement functions TSATL,TSATI are used below so that Lowe's
c         polyomial for liquid is used if T gt TMELT, or for ice if 
c         T lt TMELT. Also impose range of validity for Lowe's polys.
c
       TSATL(T) = MIN (100., MAX (T-TMELT, 0.))
       TSATI(T) = MAX (-55., MIN (T-TMELT, 0.))
c
c----------
c
c         Statement function ESAT is svp in N/m**2, with T in deg K. 
c         (100 * Lowe's poly since 1 mb = 100 N/m**2.)
c
      ESAT (T) = 
     *  100.*(
     *    CVMGT (ASAT0, BSAT0, T.GE.TMELT)
C
     *    + TSATL(T)*(ASAT1 + TSATL(T)*(ASAT2 + TSATL(T)*(ASAT3
     *    + TSATL(T)*(ASAT4 + TSATL(T)*(ASAT5 + TSATL(T)* ASAT6)))))
C
     *    + TSATI(T)*(BSAT1 + TSATI(T)*(BSAT2 + TSATI(T)*(BSAT3
     *    + TSATI(T)*(BSAT4 + TSATI(T)*(BSAT5 + TSATI(T)* BSAT6)))))
     *  )
c
c----------
c
c         Statement function DSAT is d(sat. spec. humidity)/dt, with T 
c         in deg K, and neglecting q in denominator of the q(vp) reln.
c         (100 * Lowe's poly since 1 mb = 100 N/m**2.)
c
      DSAT (T, P) = 
     *  100.*.622*(
     *    CVMGT (CSAT0, DSAT0, T.GE.TMELT)
C
     *    + TSATL(T)*(CSAT1 + TSATL(T)*(CSAT2 + TSATL(T)*(CSAT3
     *    + TSATL(T)*(CSAT4 + TSATL(T)*(CSAT5 + TSATL(T)* CSAT6)))))
C
     *    + TSATI(T)*(DSAT1 + TSATI(T)*(DSAT2 + TSATI(T)*(DSAT3
     *    + TSATI(T)*(DSAT4 + TSATI(T)*(DSAT5 + TSATI(T)* DSAT6)))))
     *  ) / P
c
c----------
c
c        Statement function QSAT is saturation specific humidity,
c        with svp E and ambient pressure P in N/m**2. Impose an upper
c        limit of 1 to avoid spurious values for very high svp
c        and/or small P.
c
       QSAT (E, P) = .622*E / MAX ( P-(1.-.622)*E, .622*E )

       QSTBLF (T, P) = .622 / MAX ( P/ESAT(T) - (1.-.622), .622 )
c
c----------