* PACKAGE DADMN  !" 力学 メインパート
*"                  ver 1.00                   92/06/12 numaguti
*"                  SRC:=$GC/model/agcm5/src/dynamics/dadmn.F
*"                  DOC:=$GC/model/agcm5/doc/dynamics.tex
**********************************************************************
      SUBROUTINE DYNMCS                      !"  力学過程制御
     M         ( GAU   , GAV   , GAT   , GAPS  , GAQ   ,
     M           GBU   , GBV   , GBT   , GBPS  , GBQ   ,
     M           GAVOR , GADIV , GBVOR , GBDIV ,
     M           GTUA  , GTVA  , GTH   , GTR   ,
     I           GTVOR , GTDIV , GTT   ,
     I           IT    , DELT  , OADVNC,
     C           ALON  , DLON  , ALAT  , DLAT  ,
     C           SIG   , SIGM  , DSIG  , DSIGM           )
*
*"              semi-implicit 時間積分
*"              implicit diffusion
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZHDIM)                      !" 文字列文字数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zhdim.F"                    !" 文字列文字数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [MODIFY]
      REAL       GAU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ
      REAL       GAV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ
      REAL       GAT   ( IDIM, JDIM, KMAX )  !" 温度  Ｔ
      REAL       GAPS  ( IDIM, JDIM       )  !" 地表気圧
      REAL       GAQ   ( IDIM, JDIM, KMAX )  !" 比湿  ｑ
      REAL       GAVOR ( IDIM, JDIM, KMAX )  !" 渦度  ζ
      REAL       GADIV ( IDIM, JDIM, KMAX )  !" 発散  Ｄ
*"       : 格子点データ(t) <GENGD> 格子点データ(t+Δt)
*
      REAL       GBU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ
      REAL       GBV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ
      REAL       GBT   ( IDIM, JDIM, KMAX )  !" 温度  Ｔ
      REAL       GBPS  ( IDIM, JDIM       )  !" 地表気圧
      REAL       GBQ   ( IDIM, JDIM, KMAX )  !" 比湿  ｑ
      REAL       GBVOR ( IDIM, JDIM, KMAX )  !" 渦度  ζ
      REAL       GBDIV ( IDIM, JDIM, KMAX )  !" 発散  Ｄ
*"       : 格子点データ(t-Δt) <MOVGD>  格子点データ(t)
*
      REAL       GTUA  ( IDIM, JDIM, KMAX )  !" 東西運動量変化項ＵＡ
      REAL       GTVA  ( IDIM, JDIM, KMAX )  !" 南北運動量変化項ＶＡ
      REAL       GTH   ( IDIM, JDIM, KMAX )  !" 温度時間変化項  Ｈ
      REAL       GTR   ( IDIM, JDIM, KMAX )  !" 比湿時間変化項  Ｒ

      REAL       GTVOR ( IDIM, JDIM, KMAX )  !" 渦度  ζ 変化
      REAL       GTDIV ( IDIM, JDIM, KMAX )  !" 発散  Ｄ 変化
      REAL       GTT   ( IDIM, JDIM, KMAX )  !" 温度  Ｔ 変化
      
*"       :  格子点変化項を足し込む
*
      INTEGER    IT                          !" 通し時間
      REAL       DELT                        !" 時間刻みΔt
      LOGICAL    OADVNC                      !" 時刻が進行するか否か
*"       :  時刻変数
*
*   [INPUT]
      REAL       ALON  ( IDIM )              !" 経度
      REAL       DLON  ( IDIM )              !" 経度荷重
      REAL       ALAT  ( JDIM )              !" 緯度
      REAL       DLAT  ( JDIM )              !" 緯度荷重
      REAL       SIG   ( KMAX )              !" σレベル(整数)
      REAL       SIGM  (KMAX+1)              !" σレベル(半整数)
      REAL       DSIG  ( KMAX )              !" Δσ(整数)
      REAL       DSIGM (KMAX+1)              !" Δσ(半整数)
*"        : 座標格子
*
*   [INTERNAL SAVE]
      REAL       DMLU  ( NMDIM, KMAX, KMAX ) !" ＳＩ行列のＬＵ分解
      INTEGER    KPIV  ( NMDIM, KMAX       ) !" ＳＩ行列のピボット
      REAL       DELTSV                      !" Δt のセーブ値
      DATA       DELTSV / 0.0  /
      SAVE       DMLU, KPIV, DELTSV
*
*   [INTERNAL ONCE]
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       EDEL  ( NMDIM  )            !" ζ，Ｄ→Ｕ，Ｖ
*
      REAL       CORIOL ( IDIM, JDIM )       !" コリオリ係数 ｆ
      REAL       UVFACT ( IDIM, JDIM )       !" u→U のファクター
*
      REAL       DIFV  ( NMDIM  )            !" 運動量水平拡散係数
      REAL       DIFT  ( NMDIM  )            !" 熱，水水平拡散係数
*
      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       WPHIS ( NMDIM )             !" 地表ジオポテンシャル
*
      SAVE       NMO   , FLAPLA, EDEL  ,
     &           CORIOL, UVFACT, DIFV  , DIFT  ,
     &           SALPHA, SBETA , SKAPPA,
     &           ATF   , BTF   , TBAR  , TBARM ,
     &           GG    , TMTX  , WHMTX , GCMTX ,
     &           WPHIS
*
      LOGICAL    OSETC , OSETV
      DATA       OSETC , OSETV  / .FALSE., .FALSE.  /
      SAVE       OSETC , OSETV
*
*   [INTERNAL PARAM]
      LOGICAL    ODIASP
      DATA       ODIASP         / .FALSE. /
      SAVE       ODIASP
*
*   [INTERNAL WORK]
      REAL       GAPIX ( IDIM, JDIM       )  !" πの東西微分
      REAL       GAPIY ( IDIM, JDIM       )  !" πの南北微分
*        : 気圧傾度
*
      REAL       GTKE  ( IDIM, JDIM, KMAX )  !" 運動エネルギー項ＫＥ
      REAL       GTUT  ( IDIM, JDIM, KMAX )  !" 温度東西移流項  ＵＴ
      REAL       GTVT  ( IDIM, JDIM, KMAX )  !" 温度南北移流項  ＶＴ
      REAL       GTPI  ( IDIM, JDIM       )  !" 気圧時間変化項  Ｚ
      REAL       GTUQ  ( IDIM, JDIM, KMAX )  !" 比湿東西移流項  Ｕｑ
      REAL       GTVQ  ( IDIM, JDIM, KMAX )  !" 比湿南北移流項  Ｖｑ
      REAL       GSIGD ( IDIM, JDIM, KMAX+1 )!" σの時間変化
      REAL       GPID  ( IDIM, JDIM       )  !" πの時間変化
*"       : 時間変化項
*
      REAL       WDVOR ( NMDIM , KMAX     )  !" 渦度  ζ
      REAL       WDDIV ( NMDIM , KMAX     )  !" 発散  Ｄ
      REAL       WDT   ( NMDIM , KMAX     )  !" 温度  Ｔ
      REAL       WDPI  ( NMDIM            )  !" log Ps＝π
      REAL       WDQ   ( NMDIM , KMAX     )  !" 比湿  ｑ
*"       : スペクトルデータ
*
      REAL       WTVOR ( NMDIM, KMAX      )  !" 渦度  ζ
      REAL       WTDIV ( NMDIM, KMAX      )  !" 発散  Ｄ
      REAL       WTT   ( NMDIM, KMAX      )  !" 温度  Ｔ
      REAL       WTPI  ( NMDIM            )  !" logPs ＝π
      REAL       WTQ   ( NMDIM, KMAX      )  !" 比湿  ｑ
      REAL       WTVDIF( NMDIM, KMAX      )  !" 渦度  ζ
      REAL       WTDDIF( NMDIM, KMAX      )  !" 発散  Ｄ

      REAL       WTVORW ( NMDIM, KMAX      )  !" 渦度  ζ
      REAL       WTDIVW ( NMDIM, KMAX      )  !" 発散  Ｄ
      REAL       WTTW   ( NMDIM, KMAX      )  !" 温度  Ｔ
      
*"       : スペクトル時間変化項
*
*
*"  << SETCON : 定数設定 >>
*
*   [ONCE]
      IF ( .NOT. OSETC ) THEN
         OSETC = .TRUE.
         CALL    DSETC                 !" スペクトル定数設定
     O         ( NMO   , FLAPLA, EDEL  ,
     O           CORIOL, UVFACT, DIFV  , DIFT  ,
     C           ALAT  , DLAT                   )
      ENDIF
*   [ONCE]
      IF ( .NOT. OSETV ) THEN
         OSETV = .TRUE.
         CALL    DSETV                 !" 鉛直差分定数設定
     O         ( SALPHA, SBETA , SKAPPA,
     O           ATF   , BTF   , TBAR  , TBARM ,
     O           GG    , TMTX  , WHMTX , GCMTX ,
     C           SIG   , SIGM  , DSIG            )
*
         CALL    DHSTRG                !" 出力設定
      ENDIF
*
      IF ( DELT .NE. DELTSV ) THEN
         CALL DMTRIX                   !" セミインプリシット用行列
     O         ( DMLU  , KPIV  ,
     I           DELT  ,
     C           WHMTX , GCMTX ,
     C           NMO   , ER    , DIFV  , DIFT  )
         DELTSV = DELT
      ENDIF
*
      CALL SETPHS                      !" 地表ジオポテンシャル設定
     O         ( WPHIS ,
     I           IT      )
*
*" << TNDNCY : 時間変化項 >>
*
      CALL PIXPIY                      !" 地表気圧傾度
     O         ( GAPIX , GAPIY ,
     I           GAPS            )
*
      CALL GRDDYN                      !"  非線型力学項
     M         ( GTUA  , GTVA  , GTH   , GTR   ,
     O           GTKE  , GTUT  , GTVT  , GTPI  , GTUQ  , GTVQ  ,
     O           GSIGD , GPID  ,
     I           GAU   , GAV   , GAT   , GAQ   ,
     I           GAVOR , GADIV , GAPIX , GAPIY ,
     C           CORIOL, UVFACT, SIGM  , DSIG  ,
     C           SALPHA, SBETA , SKAPPA,
     C           ATF   , BTF   , TBAR  , TBARM                   )
*
      CALL TENG2W                      !" 時間変化項をスペクトルに
     O         ( WTVOR , WTDIV , WTT   , WTPI  , WTQ   ,
     I           GTUA  , GTVA  , GTKE  , GTUT  , GTVT  , GTH   ,
     I           GTPI  , GTUQ  , GTVQ  , GTR   ,
     C           FLAPLA                                          )
*
      CALL GD2WD                       !" データをスペクトルに
     O         ( WDVOR , WDDIV , WDT   , WDPI  , WDQ   ,
     I           GBVOR , GBDIV , GBT   , GBPS  , GBQ    )
*

*" << GTVOR, GTDIV, GTT をスペクトルデータへ >>

         CALL DSMG2S
     O         ( WTVORW ,
     I           GTVOR  )

         CALL DSMG2S
     O         ( WTDIVW ,
     I           GTDIV  )

         CALL DSMG2S
     O         ( WTTW ,
     I           GTT  )

*" << WTVOR, WTDIV, WTT をへの追加 >>

      DO K = 1 , KMAX
         DO NM = 1 , NMDIM
            WTVOR( NM,K ) = WTVOR( NM,K ) + WTVORW( NM,K )
*
            WTDIV( NM,K ) = WTDIV( NM,K ) + WTDIVW( NM,K )
*
            WTT( NM,K )   = WTT( NM,K )   + WTTW( NM,K )
*
         END DO
      END DO
*
      
*" << INTGR : 時間積分 >>
*
      CALL TINTGR                      !" スペクトル時間積分
     M         ( WDVOR , WDDIV , WDT   , WDPI  , WDQ   ,
     I           WTVOR , WTDIV , WTT   , WTPI  , WTQ   ,
     I           DMLU  , KPIV  , DELT  ,
     C           DIFV  , DIFT  , FLAPLA, WPHIS ,
     C           GG    , TMTX  , DSIG                   )
*
      IF ( OADVNC ) THEN
         CALL MOVGD                    !" 格子点値の時刻付け変え
     O         ( GBU   , GBV   , GBT   , GBPS  , GBQ   ,
     O           GBVOR , GBDIV ,
     I           GAU   , GAV   , GAT   , GAPS  , GAQ   ,
     I           GAVOR , GADIV                          )
      ENDIF
*
      CALL GENGD                       !" 新格子点値の生成
     O         ( GAU   , GAV   , GAT   , GAPS  , GAQ   ,
     O           GAVOR , GADIV ,
     I           WDVOR , WDDIV , WDT   , WDPI  , WDQ   ,
     C           EDEL  , UVFACT                          )
*
*" << CORRCT : 補正 >>
*
      CALL TENDIF                      !" 水平拡散による時間変化
     O         ( WTVDIF, WTDDIF,
     I           WDVOR , WDDIV ,
     I           DIFV           )
*
      CALL CORDIF                      !" 水平拡散摩擦熱補正
     M         ( GAT   ,
     I           WTVDIF, WTDDIF,
     I           GAU   , GAV   , DELT  ,
     C           EDEL  , UVFACT         )
*
*" << OUTPUT : 出力 >>
*
      CALL HISTIN ( GSIGD , 'SIGDOT')  !" 出力領域へ
      CALL HISTIN ( GPID  , 'PIDOT' )
*
      CALL DDIAG                       !" 診断
     I         ( GAU   , GAV   , GAT   , GAPS  , GAQ   ,
     I           GAVOR , GADIV                          )
*
      IF ( ODIASP ) THEN
        CALL SPCTLD                    !" スペクトルでの診断
     I         ( WDVOR , 'VOR' , KMAX  , NMO   )
        CALL SPCTLD
     I         ( WDDIV , 'DIV' , KMAX  , NMO   )
        CALL SPCTLD
     I         ( WDT   , 'TMP' , KMAX  , NMO   )
        CALL SPCTLD
     I         ( WDQ   , 'Q'   , KMAX  , NMO   )
        CALL SPCTLD
     I         ( WDPI  , 'PI'  , 1     , NMO   )
      ENDIF
*
      RETURN
*======================================================================
      ENTRY DSTRT                      !" 初期データの再生成
     M         ( GAU   , GAV   , GAT   , GAPS  , GAQ   ,
     O           GAVOR , GADIV ,
     C           ALAT  , DLAT                           )
*
*   [MODIFY]
*     REAL       GAU   ( IDIM, JDIM, KMAX )  !" 西風  ｕ
*     REAL       GAV   ( IDIM, JDIM, KMAX )  !" 南風  ｖ
*     REAL       GAT   ( IDIM, JDIM, KMAX )  !" 温度  Ｔ
*     REAL       GAPS  ( IDIM, JDIM       )  !" 地表気圧
*     REAL       GAQ   ( IDIM, JDIM, KMAX )  !" 比湿  ｑ
*
*   [OUTNPUT]
*     REAL       GAVOR ( IDIM, JDIM, KMAX )  !" 渦度  ζ
*     REAL       GADIV ( IDIM, JDIM, KMAX )  !" 発散  Ｄ
*
*   [INPUT]
*     REAL       ALAT  ( JDIM )              !" 緯度
*     REAL       DLAT  ( JDIM )              !" 緯度荷重
*
*   [ONCE]
      IF ( .NOT. OSETC ) THEN
         OSETC = .TRUE.
         CALL    DSETC                 !" スペクトル定数設定
     O         ( NMO   , FLAPLA, EDEL  ,
     O           CORIOL, UVFACT, DIFV  , DIFT  ,
     C           ALAT  , DLAT                   )
      ENDIF
*
      CALL STRTV                       !" データ再生成
     M         ( GAU   , GAV   , GAT   , GAPS  , GAQ   ,
     O           GAVOR , GADIV ,
     C           EDEL  , UVFACT,
     W           WDVOR , WDDIV , WDT   , WDPI  , WDQ    )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE DHSTRG                !" 力学関係出力の登録
*
*   [INTERNAL PARAM]
      INTEGER    ISTYPL
      DATA       ISTYPL / 1 /
*
*   [INTERNAL WORK]
      REAL       VMISS
*
      CALL GZDBGT ( 'MISS' , VMISS )   !" 欠損値
*
      CALL HISTRG                      !" 出力登録
     I       ( 'SIGDOT', 'sigma-vertical vel. ' ,'1/s   ', 'AMLEV',
     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
     I         '   '   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      CALL HISTRG
     I       ( 'PIDOT ', 'pi-tendency         ' ,'1/s   ', 'ASFC',
     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
     I         '   '   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      RETURN
      END
