* PACKAGE  GTHYDRO !" T -> z hydrostatic
***********************************************************************
      PROGRAM GTHYDR
*
#ifdef SYS_IBMS
      INCLUDE    (GTSINC)
      INCLUDE    (GZSIZE)
#else
#include         "gtsinc.F"
#include         "gzsize.F"
#endif
      COMMON     /GMWORK/ MWORK
      REAL       MWORK  ( IJKDIM )
*
      CHARACTER  HHEADT( NDC )*(NCC)
      REAL       GDATAT( IJKDIM )
      CHARACTER  HHEADQ( NDC )*(NCC)
      REAL       GDATAQ( IJKDIM )
      CHARACTER  HHEADZ( NDC )*(NCC)
      REAL       GDATAZ( IJKDIM )
      CHARACTER  HITEZ  *(NCC)
*
      DATA       IFILT / 50 /
      DATA       IFILQ / 51 /
      DATA       IFILZ / 52 /
      DATA       JFILE / 60 /
*
      DATA       GRAV  /  9.8 /
*
      CHARACTER  T      *(NFILN)
      DATA       T      / 'T' /
      CHARACTER  Q      *(NFILN)
      DATA       Q      / 'q' /
      CHARACTER  ZS     *(NFILN)
      DATA       ZS     / 'ZS' /
      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  / 'Z' /
      DATA       TITLE / 'G.P. Height' /
      DATA       UNIT  / 'm' /
      DATA       DSET, EDIT, ETTL / 3*' '/
      LOGICAL    GRESET
      DATA       GRESET / .TRUE. /
      LOGICAL    HELP
      DATA       HELP   / .FALSE. /
*
      EXTERNAL   CALGPH
*     
      NAMELIST  /OPTION/ T   , Q   , ZS   , 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 ( HHEADQ, IJKDIM )
      CALL GTSIZE ( HHEADZ, IJKDIM )
      CALL GMSIZE ( IJKDIM  )
*
      CALL GURNTF ( T,  OUT  , '$GTTMPDIR/gtool.in' )
      CALL GURNTF ( Q,  OUT  , '$GTTMPDIR/gtool.in' )
      CALL GURNTF ( ZS, OUT  , '$GTTMPDIR/gtool.in' )
*
      CALL GFROPN ( IFILET, T    )
      CALL GFROPN ( IFILEQ, Q    )
      CALL GFROPN ( IFILEZ, ZS   )
      CALL GFOOPN ( JFILE ,  OUT , APND )
*
      CALL GUNENV( OUT,'.',.FALSE. )
      IL=LENC(OUT)
      WRITE (6,*) 'output='//OUT(1:IL)
*
      CALL   GFREAD
     O        ( HHEADZ, GDATAZ, IEODZ ,
     I          IFILEZ, 1               )
*
      IF ( IEODZ .NE. 0 ) THEN
         CALL MSGDMP('E','GTHYDRO','ZS FILE NOT FOUND' )
      ENDIF
*
      CALL GHCGET( HHEADZ, 'ITEM', HITEZ )
      IF ( HITEZ .EQ. 'GPHIS' ) THEN
         CALL GMFFCT
     I        ( HHEADZ, GDATAZ, 1./GRAV ,
     I          '  '  , '  '              )
      ENDIF
*
 1100 CONTINUE
         CALL   GFREAD
     O        ( HHEADT, GDATAT, IEODT ,
     I          IFILET, 1               )
         CALL   GFREAD
     O        ( HHEADQ, GDATAQ, IEODQ ,
     I          IFILEQ, 1               )
*
         IF ( MAX(IEODT,IEODQ) .EQ.0 ) THEN
            CALL GMCAL3
     I         ( CALGPH,
     M           HHEADT, GDATAT,
     I           HHEADQ, GDATAQ,
     I           HHEADZ, GDATAZ,
     I           EDIT  , ETTL  )
*
            IF ( ITEM .NE. ' ' ) THEN
               CALL GHCSET( HHEADT, 'ITEM', ITEM )
            ENDIF
            IF ( UNIT .NE. ' ' ) THEN
               CALL GHCSET( HHEADT, 'UNIT', UNIT )
            ENDIF
            IF ( TITLE .NE. ' ' ) THEN
               CALL GHCSTS( HHEADT, 'TITL', TITLE )
            ENDIF
            IF ( DSET .NE. ' ' ) THEN
               CALL GHCSET( HHEADT, 'DSET', DSET )
            ENDIF
            IF ( GRESET ) THEN
               CALL GHRSGP( HHEADT )
            ENDIF
*
            CALL  GFWRIT
     I                 ( HHEADT, GDATAT,
     I                   JFILE , 1     , 0       )
*
      GOTO 1100
         ENDIF
*
      STOP
      END
********************************************************************
      SUBROUTINE CALGPH
     I         ( HT    , T     ,
     I           HQ    , Q     ,
     I           HZS   , ZS    ,
     O           HZGP  , ZGP   ,
     D           IMAX  , JMAX  , KMAX ,
     D           IMAX2 , JMAX2 , KMAX2,
     D           IMAX3 , JMAX3 , KMAX3  )
*
      CHARACTER  HT  ( * )*(*)
      REAL       T   ( IMAX, JMAX, KMAX )
      CHARACTER  HQ  ( * )*(*)
      REAL       Q   ( IMAX2, JMAX2, KMAX3 )
      CHARACTER  HZS ( * )*(*)
      REAL       ZS  ( IMAX3, JMAX3 )
      CHARACTER  HZGP( * )*(*)
      REAL       ZGP ( IMAX, JMAX, KMAX )
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)
#else
#include         "gzsize.F"
#endif
      CHARACTER  HHS ( NDC )*(NCC)
      PARAMETER  ( KMAXD = 100 )
      REAL       SIG   ( KMAXD )
      REAL       DSIG  ( KMAXD )
*
      IF ( KMAX .GT. KMAXD ) THEN
         CALL MSGDMP( 'E','CALGPH','WORK AREA TOO SMALL' )
      ENDIF
*
      CALL GTSIZE ( HHS  ,  KMAXD   )
      CALL GUQAXV
     I            ( HT    , 3     , 'LOC',
     O              HHS   , SIG   , IEOD  )
      CALL GUQAXV
     I            ( HT    , 3     , 'WGT',
     O              HHS   , DSIG  , IEOD  )
*
      CALL HYDRO
     O         ( ZGP   ,
     I           T     , Q     , ZS    ,
     I           SIG   , DSIG  , IMAX*JMAX, KMAX  )
*
      RETURN
      END
**********************************************************************
      SUBROUTINE HYDRO
     O         ( Z     ,
     I           T     , Q     , ZS    ,
     I           SIG   , DSIG  , IJMAX , KMAX  )
*
      REAL       Z     ( IJMAX, KMAX+1 )
      REAL       T     ( IJMAX, KMAX   )
      REAL       Q     ( IJMAX, KMAX   )
      REAL       ZS    ( IJMAX         )
      REAL       SIG   ( KMAX )
      REAL       DSIG  ( KMAX )
*
      DATA       CP   / 1004.6 /  
      DATA       RAIR / 287.04 /
      DATA       GRAV / 9.8 /
      DATA       RVAP / 461. /
*
      PARAMETER  ( KMAXD = 100 )
      REAL       SALPHA ( KMAXD )
      REAL       SBETA  ( KMAXD )
*
      IF ( KMAX .GT. KMAXD ) THEN
         CALL MSGDMP( 'E','HYDRO','WORK AREA TOO SMALL' )
      ENDIF
*
      AKAPPA = RAIR / CP
      EPSV   = RAIR / RVAP
      EPSVT  = 1.0/EPSV - 1.0
*
      DO 100 K = 1, KMAX
         SALPHA( K ) = -1. +    ( 1.+ DSIG(K)/SIG(K)/2.     )**AKAPPA
         SBETA ( K ) =  1. - MAX( 1.- DSIG(K)/SIG(K)/2., 0. )**AKAPPA
  100 CONTINUE
*
      DO 1100 IJ = 1, IJMAX
           Z ( IJ,1 ) = ZS( IJ )
     &                + CP/GRAV * SALPHA( 1   ) 
     &                  * T( IJ,1 ) * ( 1.+EPSVT* Q( IJ,1 ) )
 1100 CONTINUE 
*
      DO 2300 K = 2, KMAX
         DO 2310 IJ = 1, IJMAX
            Z( IJ,K ) = Z( IJ,K-1 )
     &                + CP/GRAV * SALPHA( K   )
     &                  * T( IJ,K   ) * ( 1.+EPSVT* Q( IJ,K   ) )
     &                + CP/GRAV * SBETA ( K-1 )
     &                  * T( IJ,K-1 ) * ( 1.+EPSVT* Q( IJ,K-1 ) )
 2310    CONTINUE
 2300 CONTINUE
*
      RETURN
      END
