* PACKAGE  GTQSAT !" Qsat
***********************************************************************
      PROGRAM GTQSAT
*
#ifdef SYS_IBMS
      INCLUDE    (GTSINC)
      INCLUDE    (GZSIZE)
#else
#include         "gtsinc.F"
#include         "gzsize.F"
#endif
      CHARACTER  HHEADT( NDC )*(NCC)
      REAL       GDATAT( IJKDIM )
      CHARACTER  HHEADP( NDC )*(NCC)
      REAL       GDATAP( IJKDIM )
      CHARACTER  HHEADQ( NDC )*(NCC)
      REAL       GDATAQ( IJKDIM )
*
      DATA       IFILT / 51 /
      DATA       IFILP / 52 /
      DATA       JFILE / 60 /
*
      CHARACTER  T      *(NFILN)
      DATA       T      / 'T' /
      CHARACTER  PS     *(NFILN)
      DATA       PS     / 'Ps' /
      CHARACTER  OUT    *(NFILN)
      DATA       OUT    / '$GTTMPDIR/gtool.out' /
      LOGICAL    APND
      DATA       APND   / .FALSE. /
      CHARACTER  ITEM   *(NCC)
      CHARACTER  UNIT   *(NCC)
      CHARACTER  TITLE  *(NCC*2)
      CHARACTER  DSET   *(NCC)
      CHARACTER  EDIT   *(NCC)
      CHARACTER  ETTL   *(NCC)
      DATA       ITEM   / 'QSAT' /
      DATA       UNIT   / ' ' /
      DATA       TITLE  / 'saturation humidity' /
      DATA       DSET, EDIT, ETTL /3*' '/
      LOGICAL    GRESET
      DATA       GRESET / .FALSE. /
      LOGICAL    HELP
      DATA       HELP   / .FALSE. /
*
      EXTERNAL   QSATUM
*     
      NAMELIST  /OPTION/ T, PS, OUT, APND,
     &                   ITEM, UNIT, TITLE, DSET, EDIT, ETTL, GRESET,
     &                   HELP
*
      CALL OPTARG ( 91, 'OPTION', 'HFILE', NOPT, NFILE )
      READ (91,OPTION,IOSTAT=IOS)
      CLOSE(91)
      IF ( IOS.NE.0 .OR. HELP ) THEN
         WRITE(6,OPTION)
         STOP
      ENDIF
*
      CALL GTOPEN
      CALL GTSIZE ( HHEADT, IJKDIM )
      CALL GTSIZE ( HHEADP, IJKDIM )
      CALL GMSIZE ( IJKDIM  )
*
      CALL GURNTF ( T , OUT  , '$GTTMPDIR/gtool.in' )
      CALL GURNTF ( PS, OUT  , '$GTTMPDIR/gtool.in' )
*
      CALL GFROPN ( IFILET, T     )
      CALL GFROPN ( IFILEP, PS    )
      CALL GFOOPN ( JFILE ,  OUT , APND )
*
      CALL GUNENV( OUT,'.',.FALSE. )
      IL=LENC(OUT)
      WRITE (6,*) 'output='//OUT(1:IL)
*
 1100 CONTINUE
         CALL   GFREAD
     O        ( HHEADT, GDATAT, IEODT  ,
     I          IFILET , 1               )
         CALL   GFREAD
     O        ( HHEADP, GDATAP, IEODP ,
     I          IFILEP, 1               )
*
         IF ( MAX(IEODT,IEODP) .EQ.0 ) THEN
            CALL GPCAL2
     I         ( QSATUM,
     I           HHEADT, GDATAT,
     I           HHEADP, GDATAP,
     I           EDIT  , ETTL  ,
     O           HHEADQ, GDATAQ  )
*
            IF ( ITEM .NE. ' ' ) THEN
               CALL GHCSET( HHEADQ, 'ITEM', ITEM )
            ENDIF
            IF ( UNIT .NE. ' ' ) THEN
               CALL GHCSET( HHEADQ, 'UNIT', UNIT )
            ENDIF
            IF ( TITLE .NE. ' ' ) THEN
               CALL GHCSTS( HHEADQ, 'TITL', TITLE )
            ENDIF
            IF ( DSET .NE. ' ' ) THEN
               CALL GHCSET( HHEADQ, 'DSET', DSET )
            ENDIF
            IF ( GRESET ) THEN
               CALL GHRSGP( HHEADQ  )
            ENDIF
*
            CALL  GFWRIT
     I                 ( HHEADQ, GDATAQ,
     I                   JFILE , 1     , 0       )
*
      GOTO 1100
         ENDIF
*
      STOP
      END
********************************************************************
      SUBROUTINE QSATUM
     I         ( HHEADT, TT    ,
     I           HHEADP, PS    ,
     O           HHEADQ, QQS   ,
     D           IXDIM , IYDIM , IZDIM ,
     D           IXDIM2, IYDIM2, IZDIM2 )
*
      CHARACTER  HHEADT ( * ) *(*)
      CHARACTER  HHEADP ( * ) *(*)
      CHARACTER  HHEADQ ( * ) *(*)
      REAL       TT ( IXDIM , IYDIM , IZDIM  )
      REAL       PS ( IXDIM2, IYDIM2         )
      REAL       QQS( IXDIM , IYDIM , IZDIM  )
*
      PARAMETER  ( NCC=16, NDC=64 )
      PARAMETER  ( KDIMD=20 )
      CHARACTER  HHEADS ( NDC ) *(NCC)
      REAL       SIG ( KDIMD )
*
* PACKAGE ZQSAT     !" Clausis-Clapeyron $@$N<0(J($@J84X?t(J)
*
*"    QSAT:  saturation water vapour mixing ratio
*"    DQSAT: d(QSAT)/d(T)
*
      REAL       FQSAT
      REAL       FDQSAT
      REAL       T, P, QS
*
      FQSAT ( T,P )    = EPSV * ES0
     &                   * EXP( EL/RVAP *( 1./273. - 1./T ) )
     &                   / P
      FDQSAT( T,QS )   = EL * QS / ( RVAP * T*T )
*
      EL   = 2.5E6
      RVAP = 461.
      EPSV = 0.622
      ES0  = 6.11
      CALL GTSIZE ( HHEADS, KDIMD )
*
      CALL GUQAXV
     I         ( HHEADT, 3     , 'LOC' ,
     O           HHEADS, SIG   , IEOD   )
*
      DO 1100 IZ = 1, IZDIM
         DO 1110 IY = 1, IYDIM
            DO 1120 IX = 1, IXDIM
               QQS ( IX,IY,IZ ) 
     &            = FQSAT( TT(IX,IY,IZ), PS(IX,IY)*SIG(IZ) )
 1120       CONTINUE
 1110    CONTINUE
 1100 CONTINUE
*
      CALL GHCSET ( HHEADQ, 'ITEM', 'QSAT' )
      CALL GHCSTS ( HHEADQ, 'TITL', 'saturation humidity' )
      CALL GHCSET ( HHEADQ, 'UNIT', ' '  )
      CALL GHRSGP ( HHEADQ  )
*
      RETURN
      END
