* PACKAGE DSETV !" 力学 定数設定(鉛直差分)
*
*" 履歴 ??/??/?? 沼口敦            
*" 履歴 95/01/23 石渡正樹  TBAR0 を NAMELIST 化
************************************************************************
      SUBROUTINE DSETV       !" 力学パートの定数セット(鉛直差分)
     O         ( SALPHA, SBETA , SKAPPA,
     O           ATF   , BTF   , TBAR  , TBARM ,
     O           GG    , TMTX  , WHMTX , GCMTX ,
     C           SIG   , SIGM  , DSIG            )
*
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [OUTPUT] 
      REAL       SALPHA( KMAX )              !" 静水圧の式の係数α
      REAL       SBETA ( KMAX )              !" 静水圧の式の係数β
      REAL       SKAPPA( KMAX )              !" 温度鉛直補間の係数κ
      REAL       ATF   ( KMAX )              !" 温度鉛直補間の係数ａ
      REAL       BTF   ( KMAX )              !" 温度鉛直補間の係数ｂ
      REAL       TBAR  ( KMAX )              !" 基本温度：整数レベル
      REAL       TBARM ( KMAX )              !" 基本温度：半整数
      REAL       GG    ( KMAX )              !" セミインプリシット用
      REAL       TMTX  ( KMAX, KMAX )        !" セミインプリシット用
      REAL       WHMTX ( KMAX, KMAX )        !" セミインプリシット用
      REAL       GCMTX ( KMAX, KMAX )        !" セミインプリシット用
*"       :  鉛直差分定数
*
      REAL       SIG   ( KMAX )              !" σレベル(整数)
      REAL       SIGM  ( KMAX+1 )            !" σレベル(半整数)
      REAL       DSIG  ( KMAX )              !" Δσ
*"       :  座標格子定数
*
*   [INTERNAL PARM] 
      REAL       TBAR0                       !" 基本温度
      DATA       TBAR0 / 300.   /
*
      NAMELIST  /NMTBAR/ TBAR0
*
*   [INTERNAL WORK] 
      INTEGER    IFPAR, JFPAR
*
*
      WRITE ( 6,* ) ' DYNAMICS SETTING DATE=95/01/23'
*
      CALL   REWNML ( IFPAR , JFPAR )
      READ   ( IFPAR, NMTBAR, END=190 )
  190 WRITE  ( JFPAR, NMTBAR )
*
*"         < 1. 鉛直差分ファクター >
*
      CALL SETVDF
     O         ( SALPHA, SBETA , SKAPPA,
     O           ATF   , BTF   , TBAR  , TBARM ,
     C           SIG   , SIGM  , CP    , RAIR  ,
     I           TBAR0                           )
*
*"         < 2. semi-implicit 行列 >
*
      CALL SETMTX
     O         ( GG    , TMTX  , WHMTX , GCMTX ,
     I           SALPHA, SBETA , SKAPPA,
     I           TBAR  , TBARM,
     C           SIGM  , DSIG  , CP             )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE SETVDF      !" 鉛直差分ファクター
     O        ( SALPHA, SBETA , SKAPPA,
     O          ATF   , BTF   , TBAR  , TBARM ,
     I          SIG   , SIGM  , CP    , RAIR  ,
     I          TBAR0                            )
*
*"                 Arakawa & Suarez 鉛直差分 σ座標
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)
#else
#include        "zcdim.F"
#endif
*
*   [OUTPUT] 
      REAL       SALPHA( KMAX )              !" 静水圧の式の係数 α
      REAL       SBETA ( KMAX )              !" 静水圧の式の係数 β
      REAL       SKAPPA( KMAX )              !" 温度鉛直補間の係数κ
      REAL       ATF   ( KMAX )              !" 温度鉛直補間の係数ａ
      REAL       BTF   ( KMAX )              !" 温度鉛直補間の係数ｂ
      REAL       TBAR  ( KMAX )              !" 基本温度：整数レベル
      REAL       TBARM ( KMAX )              !" 基本温度：半整数
*
*   [INPUT] 
      REAL       SIG   ( KMAX )              !" σレベル(整数)
      REAL       SIGM  ( KMAX+1 )            !" σレベル(半整数)
      REAL       CP                          !" 大気定圧比熱
      REAL       RAIR                        !" 大気気体定数
      REAL       TBAR0                       !" 基本温度
*
*   [INTERNAL WORK] 
      REAL       AKAPPA                      !" Ｒ／Ｃｐ
      INTEGER    K
*
      AKAPPA = RAIR / CP
*
*"         < 1. 静水圧の式 >
*
      CALL HYDROS
     I   ( SIG   , SIGM  )
*
      CALL HYDROQ
     O   ( SALPHA, SBETA )
*
*"         < 2. 気圧傾度項 >
*
      DO 2100 K = 1, KMAX
         SKAPPA( K ) =  (    SIGM( K   ) * SALPHA( K )
     &                     + SIGM( K+1 ) * SBETA ( K )    )
     &                 /(    SIGM( K   ) - SIGM  ( K+1 )  )
 2100 CONTINUE
*
*"         < 3. 温度の補間 >
*
      DO 3100 K = 2, KMAX
         ATF ( K ) =   SALPHA( K )
     &              / (   1. - ( SIG( K ) / SIG( K-1 ) )**AKAPPA   )
 3100 CONTINUE
*
      DO 3200 K = 1, KMAX-1
         BTF ( K ) =   SBETA ( K )
     &              / (   ( SIG( K ) / SIG( K+1 ) )**AKAPPA - 1.   )
 3200 CONTINUE
*
*"         < 4. 基本温度場 >
*
      DO 4100 K = 1, KMAX
         TBAR  ( K ) = TBAR0
 4100 CONTINUE
*
      DO 4200 K = 2, KMAX
         TBARM ( K ) =   ATF ( K   ) * TBAR ( K   )
     &                 + BTF ( K-1 ) * TBAR ( K-1 )
 4200 CONTINUE
*
      RETURN
      END
***********************************************************************
      SUBROUTINE SETMTX          !" semi-implicit 用の行列
     O         ( GG    , TMTX  , WHMTX , GCMTX ,
     I           SALPHA, SBETA , SKAPPA,
     I           TBAR  , TBARM ,
     I           SIGM  , DSIG  , CP             )
*
*"                 Arakawa & Suarez 鉛直差分 σ座標
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)
#else
#include        "zcdim.F"
#endif
*
*   [OUTPUT] 
      REAL       GG    ( KMAX )              !" セミインプリシット用
      REAL       TMTX  ( KMAX, KMAX )        !" セミインプリシット用
      REAL       WHMTX ( KMAX, KMAX )        !" セミインプリシット用
      REAL       GCMTX ( KMAX, KMAX )        !" セミインプリシット用
*
*   [INPUT] 
      REAL       SALPHA( KMAX )              !" 静水圧の式の係数 α
      REAL       SBETA ( KMAX )              !" 静水圧の式の係数 β
      REAL       SKAPPA( KMAX )              !" 温度鉛直補間の係数κ
      REAL       TBAR  ( KMAX )              !" 基本温度：整数レベル
      REAL       TBARM ( KMAX )              !" 基本温度：半整数
      REAL       SIGM  ( KMAX+1 )            !" σレベル(半整数)
      REAL       DSIG  ( KMAX )              !" Δσ
      REAL       CP                          !" 大気定圧比熱
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           WW    , SS    ,
     &           QQ    , RR
      REAL       WW    ( KMAX, KMAX )        !" Ｗ
      REAL       SS    ( KMAX, KMAX )        !" Ｓ
      REAL       QQ    ( KMAX, KMAX )        !" Ｑ
      REAL       RR    ( KMAX, KMAX )        !" Ｒ
*
      INTEGER    K, L, M
*
*"         < 0. Ｃ = Δσ >
*
*"         < 1. Ｇ = ＣｐκＴ >
*
      DO 1100 K = 1, KMAX
         GG ( K ) = CP * SKAPPA( K ) * TBAR( K )
 1100 CONTINUE
*
*"         < 2. Ｗ : Φ = Ｗ Ｔ >
*
      CALL RESET ( WW, KMAX*KMAX )
*
      DO 2120 K = 1, KMAX
*
         DO 2100 L = 1, K
            WW( K,L ) =             CP * SALPHA( L )
 2100   CONTINUE
*
         DO 2110 L = 1, K-1
            WW( K,L ) = WW( K,L ) + CP * SBETA( L )
 2110    CONTINUE
*
 2120 CONTINUE
*
*"         < 3. Ｓ ： ｄσ／ｄｔ ＝ Ｓ Ｄ >
*
      DO 3120 K = 1, KMAX
*
         DO 3100 L = 1, KMAX
            SS( K,L ) =   SIGM( K ) * DSIG( L )
 3100    CONTINUE
*
         DO 3110 L = K, KMAX
            SS( K,L ) =   SS( K,L ) - DSIG( L )
 3110    CONTINUE
*
 3120 CONTINUE
*
*"         < 4. Ｑ ： ｄσ／ｄｔ ∂Ｔ/∂σ ＝ Ｑ Ｓ Ｄ >
*
      CALL RESET ( QQ, KMAX*KMAX )
*
      DO 4100 K = 1, KMAX
         QQ( K,K   ) = ( TBARM( K ) - TBAR ( K   ) ) / DSIG( K )
 4100 CONTINUE
*
      DO 4110 K = 1, KMAX-1
         QQ( K,K+1 ) = ( TBAR ( K ) - TBARM( K+1 ) ) / DSIG( K )
 4110 CONTINUE
*
*"         < 5. Ｒ：κＴ（∂π/∂ｔ＋ｄσ/ｄｔ／σ）
*"                  ＝ Ｒ Ｄ >
*
      CALL RESET ( RR, KMAX*KMAX )
*
      DO 5120 K = 1, KMAX
*
         DO 5100 L = K, KMAX
            RR( K,L ) = - SALPHA( K ) / DSIG( K ) * DSIG( L )
     &                    * TBAR( K )
 5100    CONTINUE
*
         DO 5110 L = K+1, KMAX
            RR( K,L ) = RR( K,L )
     &                   - SBETA ( K ) / DSIG( K ) * DSIG( L )
     &                     * TBAR( K )
 5110    CONTINUE
*
 5120 CONTINUE
*
*"         < 6. Ｑ Ｓ >
*
      CALL RESET( TMTX, KMAX*KMAX )
*
      DO 6100 K = 1, KMAX
         DO 6100 L = 1, KMAX
            DO 6100 M = 1, KMAX
               TMTX( K,L ) = TMTX( K,L ) + QQ( K,M ) * SS( M,L )
 6100 CONTINUE
*
*"         < 7. ｈ ＝ Ｑ Ｓ - Ｒ >
*
      DO 7100 K = 1, KMAX
         DO 7100 L = 1, KMAX
            TMTX( K,L ) = TMTX( K,L ) - RR( K,L )
 7100 CONTINUE
*
*"         < 8. Ｗ ｈ >
*
      CALL RESET ( WHMTX, KMAX*KMAX )
*
      DO 8100 K = 1, KMAX
         DO 8100 L = 1, KMAX
            DO 8100 M = 1, KMAX
               WHMTX( K,L ) = WHMTX( K,L ) + WW( K,M ) * TMTX( M,L )
 8100 CONTINUE
*
*"         < 9. Ｇ tＣ  >
*
      DO 9100 K = 1, KMAX
         DO 9100 L = 1, KMAX
               GCMTX( K,L ) = GG( K ) * DSIG( L )
 9100 CONTINUE
*
      RETURN
      END
