*%表題   静的エネルギー, エントロピー関連の諸量を計算
*%
*%履歴   91/05/11 佐藤正樹       C:\FORT\OORT\MKSEM.FOR
*%       91/05/25 佐藤正樹       C:\FORT\OORT\MKSEER.FOR
*%
*%
*-----------------------------------------------------------------------
      PROGRAM  MKSEER

      INTEGER  NLAT           ! 緯度配列の次元
      PARAMETER ( NLAT = 18 )
      INTEGER  NPRS           ! 圧力配列の次元
      PARAMETER ( NPRS = 20 )

      REAL  PRS(0:NPRS)      ! 圧力
      REAL  TMP(-NLAT:NLAT, 0:NPRS) ! 温度
      REAL  QHU(-NLAT:NLAT, 0:NPRS) ! 比湿
      REAL  GPZ(-NLAT:NLAT, 0:NPRS) ! ジオポテンシャル高度

      REAL  SEM(-NLAT:NLAT, 0:NPRS) ! 湿潤静的エネルギー
      REAL  SED(-NLAT:NLAT, 0:NPRS) ! 乾燥静的エネルギー
      REAL  ERM(-NLAT:NLAT, 0:NPRS) ! 湿潤エントロピー
      REAL  ERD(-NLAT:NLAT, 0:NPRS) ! 乾燥エントロピー
      REAL  PTM(-NLAT:NLAT, 0:NPRS) ! 温位
      REAL  QRE(-NLAT:NLAT, 0:NPRS) ! 相対湿度



      REAL*8  R                     ! 惑星半径
      REAL*8  OMEGA0                ! 自転角速度
      DOUBLE PRECISION  PI          ! 円周率
      PARAMETER( PI = 3.141592653589793 D 0 )



      INTEGER  ILAT           ! 緯度配列の次元, DOループ
      INTEGER  IPRS           ! 圧力配列の次元, DOループ

      REAL*8   PRS2           ! 倍精度圧力
      REAL*8   TMP2           ! 倍精度温度
      REAL*8   QHU2           ! 倍精度比湿
      REAL*8   GPZ2           ! 倍精度ジオポテンシャル高度



*     < 関数プログラム >

      REAL*8    FSEMQ         ! 湿潤静的エネルギー
      REAL*8    TFETRQ        ! エントロピー
      REAL*8    TFPTMD        ! 温位
      REAL*8    TFQREQ        ! 相対湿度



*     < 入出力パラメーター >

      INTEGER   IPR  / 10 /          ! 入力ユニット番号
      CHARACTER HIPFL*13 / 'MKSEER.PAR' / 
     &                               ! 入力パラメーターファイル名
      LOGICAL   LIPEX                ! 入力パラメーターファイルの存在

      INTEGER   IR   / 20 /          ! 入力ユニット番号
      INTEGER  IIFL                  ! 入力ファイル数
      PARAMETER( IIFL=4 )
      CHARACTER HIFL(IIFL)*16        ! 入力ファイル名

      INTEGER   OR   / 30 /          ! 出力ユニット番号
      INTEGER  IOFL                  ! 出力ファイル数
      PARAMETER( IOFL=6 )
      CHARACTER HOFL(IOFL)*16        ! 出力ファイル名

      CHARACTER  HXFORM*16           ! 実数型変数の入出力形式



*     < GTOOL >

      INTEGER  NCC
      INTEGER  NDC
      PARAMETER  ( NCC=16, NDC=64 )
*
      CHARACTER  HHEAD ( NDC )*(NCC)
*
      CHARACTER  HDSET *(NCC)
      CHARACTER  HITEM *(NCC)
      CHARACTER  HTITL *(NCC*2)
      CHARACTER  HUNIT *(NCC)
      CHARACTER  HUTIM *(NCC)
      CHARACTER  HAITM1*(NCC)
      CHARACTER  HAITM2*(NCC)
      CHARACTER  HAITM3*(NCC)
      CHARACTER  HEDIT*(NCC)
      CHARACTER  HETTL*(NCC)

      INTEGER    IFNUM
      INTEGER    ITDUR
      REAL       DMIN
      REAL       DMAX
      REAL       DIVS
      REAL       DIVL
      INTEGER    ISTYP

      INTEGER    IYEAR
      INTEGER    IMONTH
      INTEGER    IDAY
      INTEGER    IHOUR
      INTEGER    IMIN
      INTEGER    ISEC

      INTEGER  IDIM
      INTEGER  JDIM
      INTEGER  KDIM
      INTEGER  IJKDIM



*     < NAMELIST >

      NAMELIST / FLPARA /  HIFL,
     &                     HOFL

*-----------------------------------------------------------------------
*     < アンダーフローのスイッチの設定 >

      CALL UNDER0(.TRUE.)

*-----------------------------------------------------------------------
*     < 入力形式 >

      HXFORM='(5E16.6)'

*-----------------------------------------------------------------------
*     < 定数 >

      R      = 6378.D3
      OMEGA0 = 2.D0 * PI / 24.D0 / 60.D0 / 60.D0

*-----------------------------------------------------------------------
*     < 入力データファイル名の入力 >

      INQUIRE(FILE=HIPFL, EXIST=LIPEX)
      IF( .NOT. LIPEX ) THEN
        WRITE(6,*) '入力パラメーターファイルが存在しません'
        WRITE(6,*) ' ファイル名: ', HIPFL
        STOP
      END IF


      OPEN(IPR,FILE=HIPFL)
        READ(IPR,FLPARA)
      CLOSE(IPR)

      WRITE(6,*) 'ファイルを読み込みます', HIFL

*-----------------------------------------------------------------------
*     < データ入力 >

      CALL  RDGTDT
     I           ( IR,       HIFL(1),
     O             HHEAD,    PRS,
     D             IDIM,     JDIM,     KDIM,     IJKDIM  )

      IF( IDIM .NE. NPRS+1 ) THEN

         WRITE(6,*) '配列の数がおかしい', HIFL(1)
         STOP
       ENDIF



      CALL  RDGTDT
     I           ( IR,       HIFL(2),
     O             HHEAD,    TMP,
     D             IDIM,     JDIM,     KDIM,     IJKDIM  )

      IF( ( IDIM .NE. 1 ) .OR. ( JDIM .NE. (2*NLAT+1) ) 
     &  .OR.  ( KDIM .NE. NPRS+1 ) 
     &  .OR. ( IJKDIM .NE. (2*NLAT+1) * ( NPRS+1) ) ) THEN

         WRITE(6,*) '配列の数がおかしい', HIFL(2)
         STOP
       ENDIF


      CALL  RDGTDT
     I           ( IR,       HIFL(3),
     O             HHEAD,    QHU,
     D             IDIM,     JDIM,     KDIM,     IJKDIM  )

      IF( ( IDIM .NE. 1 ) .OR. ( JDIM .NE. (2*NLAT+1) ) 
     &  .OR.  ( KDIM .NE. NPRS+1 ) 
     &  .OR. ( IJKDIM .NE. (2*NLAT+1) * ( NPRS+1) ) ) THEN

         WRITE(6,*) '配列の数がおかしい', HIFL(3)
         STOP
       ENDIF


      CALL  RDGTDT
     I           ( IR,       HIFL(4),
     O             HHEAD,    GPZ,
     D             IDIM,     JDIM,     KDIM,     IJKDIM  )

      IF( ( IDIM .NE. 1 ) .OR. ( JDIM .NE. (2*NLAT+1) ) 
     &  .OR.  ( KDIM .NE. NPRS+1 ) 
     &  .OR. ( IJKDIM .NE. (2*NLAT+1) * ( NPRS+1) ) ) THEN

         WRITE(6,*) '配列の数がおかしい', HIFL(4)
         STOP
       ENDIF

*-----------------------------------------------------------------------
*     < 湿潤静的エネルギーの計算 >

      CALL TSCERT                    ! 熱力学ライブラリ

      DO 90 ILAT = -NLAT, NLAT
        SEM(ILAT, 0) =  -999.
        SED(ILAT, 0) =  -999.
        ERM(ILAT, 0) =  -999.
        ERD(ILAT, 0) =  -999.
        PTM(ILAT, 0) =  -999.
        QRE(ILAT, 0) =  -999.
 90   CONTINUE

      DO 100 IPRS = 1, NPRS
        PRS2 = DBLE( PRS(IPRS) )
        DO 200 ILAT = -NLAT, NLAT
          TMP2 = DBLE( TMP(ILAT, IPRS) )
          QHU2 = DBLE( QHU(ILAT, IPRS) )
          GPZ2 = DBLE( GPZ(ILAT, IPRS) )
          SEM(ILAT, IPRS) =  REAL( FSEMQ( PRS2, TMP2, QHU2, GPZ2 ) )
          SED(ILAT, IPRS) =  REAL( FSEMQ( PRS2, TMP2, 0.D0, GPZ2 ) )
          ERM(ILAT, IPRS) =  REAL( TFETRQ( PRS2, TMP2, QHU2 ) )
          ERD(ILAT, IPRS) =  REAL( TFETRQ( PRS2, TMP2, 0.D0 ) )
          PTM(ILAT, IPRS) =  REAL( TFPTMD( PRS2, TMP2 ) )
          QRE(ILAT, IPRS) =  REAL( TFQREQ( PRS2, TMP2, QHU2 ) )
 200  CONTINUE
 100  CONTINUE

*-----------------------------------------------------------------------
*     < パラメーターの設定 >

      HDSET  = 'Oort and Peixoto'

      IFNUM  = 1
      HUTIM  = 'SEC'
      ITDUR  = 31536000.
      HAITM1 = 'GLON1'
      HAITM2 = 'GLAT37'
      HAITM3 = 'GPRS21'

      HEDIT  = 'YR'
      HETTL  = 'YEAR 63-73'

      IYEAR  = 1963
      IMONTH = 0
      IDAY   = 0
      IHOUR  = 0
      IMIN   = 0
      ISEC   = 0

*-----------------------------------------------------------------------
*     < SEM の出力 >

      HITEM  = 'SEM'
      HTITL  = 'Moist Static Energy'
      HUNIT  = 'J/kg'

      DMIN   = 2.0E5
      DMAX   = 5.0E5
      DIVS   = 0.1E5
      DIVL   = 0.5E5
      ISTYP  = 1

      CALL GTDATA
     I           ( OR,     HOFL(1),    HXFORM, 
     I             SEM,
     I             HDSET , HITEM , IFNUM ,
     I             HTITL , HUNIT , HUTIM , ITDUR ,
     I             HEDIT,  HETTL, 
     I             HAITM1, HAITM2, HAITM3,
     I             DMIN  , DMAX  , DIVS  , DIVL  , ISTYP,
     I             IYEAR , IMONTH, IDAY  , IHOUR , IMIN  , ISEC ,
     D             1,      2*NLAT+1, NPRS+1 )

*-----------------------------------------------------------------------
*     < SED の出力 >

      HITEM  = 'SED'
      HTITL  = 'Dry Static Energy'
      HUNIT  = 'J/kg'

      DMIN   = 2.0E5
      DMAX   = 5.0E5
      DIVS   = 0.1E5
      DIVL   = 0.5E5
      ISTYP  = 1

      CALL GTDATA
     I           ( OR,     HOFL(2),    HXFORM, 
     I             SED,
     I             HDSET , HITEM , IFNUM ,
     I             HTITL , HUNIT , HUTIM , ITDUR ,
     I             HEDIT,  HETTL, 
     I             HAITM1, HAITM2, HAITM3,
     I             DMIN  , DMAX  , DIVS  , DIVL  , ISTYP,
     I             IYEAR , IMONTH, IDAY  , IHOUR , IMIN  , ISEC ,
     D             1,      2*NLAT+1, NPRS+1 )

*-----------------------------------------------------------------------
*     < ERM の出力 >

      HITEM  = 'ERM'
      HTITL  = 'Moist Entropy'
      HUNIT  = 'J/kg K'

      DMIN   = -100.
      DMAX   = 500.
      DIVS   = 20.
      DIVL   = 100.
      ISTYP  = 1

      CALL GTDATA
     I           ( OR,     HOFL(3),    HXFORM, 
     I             ERM,
     I             HDSET , HITEM , IFNUM ,
     I             HTITL , HUNIT , HUTIM , ITDUR ,
     I             HEDIT,  HETTL, 
     I             HAITM1, HAITM2, HAITM3,
     I             DMIN  , DMAX  , DIVS  , DIVL  , ISTYP,
     I             IYEAR , IMONTH, IDAY  , IHOUR , IMIN  , ISEC ,
     D             1,      2*NLAT+1, NPRS+1 )

*-----------------------------------------------------------------------
*     < ERD の出力 >

      HITEM  = 'ERD'
      HTITL  = 'Dry Entropy'
      HUNIT  = 'J/kg K'

      DMIN   = -100.
      DMAX   = 500.
      DIVS   = 20.
      DIVL   = 100.
      ISTYP  = 1

      CALL GTDATA
     I           ( OR,     HOFL(4),    HXFORM, 
     I             ERD,
     I             HDSET , HITEM , IFNUM ,
     I             HTITL , HUNIT , HUTIM , ITDUR ,
     I             HEDIT,  HETTL, 
     I             HAITM1, HAITM2, HAITM3,
     I             DMIN  , DMAX  , DIVS  , DIVL  , ISTYP,
     I             IYEAR , IMONTH, IDAY  , IHOUR , IMIN  , ISEC ,
     D             1,      2*NLAT+1, NPRS+1 )

*-----------------------------------------------------------------------
*     < PTM の出力 >

      HITEM  = 'PTM'
      HTITL  = 'Potential Temperature'
      HUNIT  = 'K'

      DMIN   = 200.
      DMAX   = 500.
      DIVS   =  5.
      DIVL   =  20.
      ISTYP  =   1

      CALL GTDATA
     I           ( OR,     HOFL(5),    HXFORM, 
     I             PTM,
     I             HDSET , HITEM , IFNUM ,
     I             HTITL , HUNIT , HUTIM , ITDUR ,
     I             HEDIT,  HETTL, 
     I             HAITM1, HAITM2, HAITM3,
     I             DMIN  , DMAX  , DIVS  , DIVL  , ISTYP,
     I             IYEAR , IMONTH, IDAY  , IHOUR , IMIN  , ISEC ,
     D             1,      2*NLAT+1, NPRS+1 )

*-----------------------------------------------------------------------
*     < QRE の出力 >

      HITEM  = 'QRE'
      HTITL  = 'Relative Humidity'
      HUNIT  = ' '

      DMIN   = 0.
      DMAX   = 1.
      DIVS   = 0.1
      DIVL   = 0.2
      ISTYP  =   1

      CALL GTDATA
     I           ( OR,     HOFL(6),    HXFORM, 
     I             QRE,
     I             HDSET , HITEM , IFNUM ,
     I             HTITL , HUNIT , HUTIM , ITDUR ,
     I             HEDIT,  HETTL, 
     I             HAITM1, HAITM2, HAITM3,
     I             DMIN  , DMAX  , DIVS  , DIVL  , ISTYP,
     I             IYEAR , IMONTH, IDAY  , IHOUR , IMIN  , ISEC ,
     D             1,      2*NLAT+1, NPRS+1 )

*-----------------------------------------------------------------------
      STOP
      END
*=======================================================================
*%表題   湿潤静的エネルギーを求める
*%
*%履歴   91/05/11  佐藤正樹     C:\FORT\OORT\MKSEM.FOR
*%       91/05/25  佐藤正樹     C:\FORT\OORT\MKSEER.FOR
*%
*-----------------------------------------------------------------------
      FUNCTION  FSEMQ
     I        ( PRS,      TMP,      QHU,      ZZZ )
*
      REAL*8    FSEMQ    ! 湿潤静的エネルギー

      REAL*8    PRS      ! 圧力
      REAL*8    TMP      ! 温度
      REAL*8    QHU      ! 比湿
      REAL*8    ZZZ      ! 高度

*     < SUBROUTINE 内で定義される変数 >

      REAL*8    TFETHQ   ! エンタルピー

*-----------------------------------------------------------------------
*     < 定数の引渡し :
*       ファイル名 TICNST, 共通ブロック名 / CMCNST / >

      INCLUDE 'TICNST'

*-----------------------------------------------------------------------
*     < 温位 >

      FSEMQ = TFETHQ( PRS,      TMP,      QHU ) + GRAV * ZZZ

*-----------------------------------------------------------------------

      RETURN
      END
*=======================================================================
*%表題   温位を求める
*%
*%履歴   91/05/11  佐藤正樹     C:\FORT\OORT\MKPTM.FOR
*%
*-----------------------------------------------------------------------
      FUNCTION TFPTMD
     I        ( PRS,      TMP )
*
      REAL*8    TFPTMD   ! 温位

      REAL*8    PRS      ! 圧力
      REAL*8    TMP      ! 温度

*     < SUBROUTINE 内で定義される変数 >

      REAL*8    TFETRQ   ! エントロピー

*-----------------------------------------------------------------------
*     < 定数の引渡し :
*       ファイル名 TICNST, 共通ブロック名 / CMCNST / >

      INCLUDE 'TICNST'

*-----------------------------------------------------------------------
*     < 温位 >

      TFPTMD = TMP0 * EXP( TFETRQ( PRS,      TMP,      0.D0 ) / CPD )

*-----------------------------------------------------------------------

      RETURN
      END
