SUBROUTINE QSAT (TT,P,Q,DQDT,LDQDT)

      _RL TT, P, Q, DQDT

      D = (P-ERFAC*QX)

      SUBROUTINE SHC2GRID( LMAX, SHC, NLON, NLAT, GRID, P )

      REAL    P((LMAX+1)*(LMAX+2)/2)     ! work space

         CALL SHC2ZONE( LMAX, SHC, XLAT1, XLAT2, HS, HN, P )

      SUBROUTINE SHC4GRID( LMAX, SHC, NLON, NLAT, GRID, P )

      REAL    P((LMAX+1)*(LMAX+2)/2)     ! work space

         CALL SHC4ZONE(LMAX,SHC,XLAT1,XLAT2,HS,HN,P)

      SUBROUTINE SHC2ZONE(LMAX,SHC,XLAT1,XLAT2,HS,HN,P)

      REAL    P((LMAX+1)*(LMAX+2)/2) ! INPUT

      CALL HELMHOLTZ(LMAX,SLAT,P)

            Q    = FACT(L)*P(K)

            Q    = FACT(L+1)*P(K)

            Q    = FACT(LMAX)*P(K)

      SUBROUTINE SHC4ZONE(LMAX,SHC,XLAT1,XLAT2,HS,HN,P)

      REAL    P((LMAX+1)*(LMAX+2)/2)

      CALL HELMHOLTZ(LMAX,SIN0,P)

            SHC(K+M) = SHC(K+M) + P(I1) * CE

            SHC(K-M) = SHC(K-M) + P(I1) * SE

            SHC(K+M) = SHC(K+M) + P(I2) * CO                    

            SHC(K-M) = SHC(K-M) + P(I2) * SO

            SHC(K+M) = SHC(K+M) + P(I) * CE

            SHC(K-M) = SHC(K-M) + P(I) * SE

      SUBROUTINE HELMHOLTZ(LMAX,S,P)

      REAL    P(1)        ! OUTPUT assoc. legendre polynomials.

         P(1) = 1.0

            P(K) = Z(M) * CM

            P(K) = S*P(K-1)/AK

               P(K) =(S*P(K-1)-P(K-2)*A(M))/AK

     &     (5.2787E-8*T-6.12293E-6)*T+3.47718E-5) *P + 

     &     *P + (((-5.3009E-4*T+1.6483E-2)*T+7.944E-2) * SR + 

      B  = (5.03217E-5*P+3.359406)*P+21582.27

      x = x*(1.0 - P/B) + (V350P+x)*P*SR/(B*(B+SR))

      SR = V350P*(1.0 - P/B)

      SBO_RHO = 1028.106331 + P/B/SR - x / (SR*(SR+x))

      _RL PLAT,D,C1,P,PR,Q,X,SR,V350P,B

      P=((1-C1)-sqrt(((1-C1)**2)-(8.84E-6*abs(DPT))))/4.42E-6

      PR = P

      P  = 0.

      Q = PR*((((-2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P+

     &     1.8741E-8))*P+(-4.2393E-8*T+1.8932E-6)*(S-35.0)+

      P = P + 0.5*PR

      x = PR*((((-2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P+

     &     1.8741E-8))*P+(-4.2393E-8*T+1.8932E-6)*(S-35.0)+

      x = PR*((((-2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P+

     &     1.8741E-8))*P+(-4.2393E-8*T+1.8932E-6)*(S-35.0)+

      P = P + 0.5*PR 

      x = PR*((((-2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P+

     &     1.8741E-8))*P+(-4.2393E-8*T+1.8932E-6)*(S-35.0)+

      P=P/10.

      _RL S,T,P,PR

      del_P  = PR - P

      del_th = del_P*sw_adtg(S,T,P)

      del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)

      del_th = del_P*sw_adtg(S,th,P+onehalf*del_P)

      del_th = del_P*sw_adtg(S,th,P+del_P)

      _RL function SW_PTMP  (S,T,P,PR)

      _RL S,T,P

     &     + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P 

     &     + (  e0 + (e1 + e2*T)*T )*P*P

      _RL function SW_ADTG  (S,T,P)