* PACKAGE DADMN  !" 力学 メインパート
*"                  ver 1.02                   93/06/23 takepiro
**********************************************************************
      SUBROUTINE DYNMCS                      !"  力学過程制御
     M         ( GAU   , GAV   , GAW   , GAT   , 
     M           GATOR , GAPOR , 
     I           IT    , DELT  , OADVNC,
     C           ALON  , DLON  , ALAT  , DLAT  ,
     C           ARAD  , DRAD                    )
*
*"      Runge-Kutta 時間積分( 水平スペクトル - 鉛直グリッド空間 )
*
*   [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, 0:KDIM )  !" 西風   ｕ
      REAL       GAV   ( IDIM, JDIM, 0:KDIM )  !" 南風   ｖ
      REAL       GAW   ( IDIM, JDIM, 0:KDIM )  !" 鉛直風 ｗ
      REAL       GAT   ( IDIM, JDIM, 0:KDIM )  !" 温度   Ｔ
      REAL       GATOR ( IDIM, JDIM, 0:KDIM )  !" トロイダル Ψ
      REAL       GAPOR ( IDIM, JDIM, 0:KDIM )  !" ポロイダル Φ
*"       : 格子点データ(t) <DYNMCS>  格子点データ(t+Δt)
*
      INTEGER    IT                          !" 通し時間
      REAL       DELT                        !" 時間刻みΔt
      LOGICAL    OADVNC                      !" 時刻が進行するか否か
*"       :  時刻変数
*
*   [INPUT]
      REAL       ALON  ( IDIM )              !" 経度
      REAL       DLON  ( IDIM )              !" 経度荷重
      REAL       ALAT  ( JDIM )              !" 緯度
      REAL       DLAT  ( JDIM )              !" 緯度荷重
      REAL       ARAD  ( 0:KDIM )            !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )            !" Δｒ(整数)
*"        : 座標格子
*
*   [INTERNAL SAVE]
*
*   [INTERNAL ONCE]
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       EDEL  ( NMDIM  )            !" ζ，Ｄ→Ｕ，Ｖ
*
      REAL       UVFACT ( IDIM, JDIM )       !" u→U のファクター
*
      REAL       RAYL                          !" レイリー数
      REAL       TAU                           !" √テイラー数
      REAL       PRND                          !" プランドル数
*
      SAVE       NMO   , FLAPLA, EDEL  ,
     &           UVFACT
*
*
*   [INTERNAL PARAM]
      LOGICAL    OSETC , OFIRST
      DATA       OSETC , OFIRST  / .FALSE., .TRUE.  /
      SAVE       OSETC , OFIRST
*
      LOGICAL    ONONLN                        !" 非線形力学スイッチ
      SAVE       ONONLN
*
      LOGICAL    ODIASP
      DATA       ODIASP         / .TRUE. /
      SAVE       ODIASP
*
*   [INTERNAL WORK]
      REAL       GTUA  ( IDIM, JDIM, 0:KDIM )  !" 東西運動量変化項ＵＡ
      REAL       GTVA  ( IDIM, JDIM, 0:KDIM )  !" 南北運動量変化項ＶＡ
      REAL       GTWA  ( IDIM, JDIM, 0:KDIM )  !" 鉛直運動量変化項ＶＡ
*
      REAL       GTUT  ( IDIM, JDIM, 0:KDIM )  !" 温度東西移流項  ＵＴ
      REAL       GTVT  ( IDIM, JDIM, 0:KDIM )  !" 温度南北移流項  ＶＴ
      REAL       GTWT  ( IDIM, JDIM, 0:KDIM )  !" 温度鉛直移流項  ＷＴ
      REAL       GTH   ( IDIM, JDIM, 0:KDIM )  !" 温度時間変化項  Ｈ
*"       : 格子点変化項を足し込む
*"       : 時間変化項
*
      REAL       WDTOR ( NMDIM , 0:KDIM     )  !" トロイダル Ψ
      REAL       WDPOR ( NMDIM , 0:KDIM     )  !" ポロイダル Φ
      REAL       WDDPOR( NMDIM , 0:KDIM     )  !" ポロイダル D_l Φ
      REAL       WDT   ( NMDIM , 0:KDIM     )  !" 温度  Ｔ
      SAVE       WDTOR , WDPOR , WDDPOR, WDT
*"       : スペクトルデータ
*
      REAL       WTTOR ( NMDIM, 0:KDIM      )  !" トロイダル Ψ
      REAL       WTPOR ( NMDIM, 0:KDIM      )  !" ポロイダル Φ
      REAL       WTT   ( NMDIM, 0:KDIM      )  !" 温度  Ｔ
*"       : スペクトル時間変化項
*
      INTEGER    ITRUNG                        !" Runge-Kutta 制御変数
      INTEGER    IORUNG                        !" Runge-Kutta 次数
      PARAMETER( IORUNG = 4 )
*
*
*"  << SETCON : 定数設定 >>
*
*   [ONCE]
      IF ( .NOT. OSETC ) THEN
         OSETC = .TRUE.
         CALL    DSETC                 !" スペクトル定数設定
     O         ( NMO   , FLAPLA, EDEL  , UVFACT, 
     C           ALAT  , DLAT                      )
      ENDIF
*
*   [ONCE]
      IF ( OFIRST ) THEN
         OFIRST = .FALSE.
         CALL    DPARM                     !" パラメター設定
     O             ( RAYL   , TAU  , PRND  , 
     O               ONONLN                   )
*
         CALL HSPARM
     I          ( 'RAYLEIGH NUMBER', RAYL   )
         CALL HSPARM
     I          ( 'TAYLOR NUMBER'  , TAU**2 )
         CALL HSPARM
     I          ( 'PRANDTL NUMBER' , PRND   )
*
         CALL    DHSTRG                    !" 出力設定
*
         CALL GD2WD             !" データをスペクトルに
     O         ( WDTOR , WDPOR , WDT  , WDDPOR ,
     I           GATOR , GAPOR , GAT  ,
     C           FLAPLA,
     C           ARAD  , DRAD                     )
*
      ENDIF
*
*" << Runge-Kutta Loop >>
*
      ITRUNG = 0
*
 8888 CONTINUE                    !" <---  Runge-Kutta Loop : ここから
        ITRUNG = ITRUNG + 1 
*
         CALL DLOPER
     M           ( WDDPOR, 
     I             WDPOR ,
     F             'POS ',
     C             FLAPLA,
     C             ARAD  , DRAD    )
*
*
*" << TNDNCY : 時間変化項 >>
*
         IF( ONONLN )THEN
            CALL GRDDYN         !"  非線型力学項
     O         ( GTUA  , GTVA  , GTWA  , 
     O           GTUT  , GTVT  , GTWT  , GTH   ,
     I           GAU   , GAV   , GAW   , GAT   ,
     C           UVFACT, FLAPLA, 
     C           ALAT  , DLAT  , ARAD  , DRAD    )
*
         ELSE
            CALL RESET( GTUA , IDIM*JDIM*KMAX )
            CALL RESET( GTVA , IDIM*JDIM*KMAX )
            CALL RESET( GTWA , IDIM*JDIM*KMAX )
            CALL RESET( GTUT , IDIM*JDIM*KMAX )
            CALL RESET( GTVT , IDIM*JDIM*KMAX )
            CALL RESET( GTWT , IDIM*JDIM*KMAX )
            CALL RESET( GTH  , IDIM*JDIM*KMAX )
         ENDIF
*
         CALL GRDFRC            !" 力学外部強制項(格子点)
     M         ( GTUA  , GTVA  , GTWA  , 
     M           GTUT  , GTVT  , GTWT  , GTH   ,
     C           ALON  , DLON  , ALAT  , DLAT  ,
     C           ARAD  , DRAD                    )
*
         CALL TENG2W            !" 時間変化項をスペクトルに
     O         ( WTTOR , WTPOR , WTT   , 
     I           GTUA  , GTVA  , GTWA  , 
     I           GTUT  , GTVT  , GTWT  , GTH , 
     C           FLAPLA,
     C           ARAD  , DRAD                  )
*
        CALL LNDYN              !" 力学 線型項
     M         ( WTTOR , WTPOR , WTT   , 
     I           WDTOR , WDPOR , WDT   , WDDPOR, 
     C           RAYL  , TAU   , PRND  , 
     C           NMO   , FLAPLA, 
     C           ARAD  , DRAD                    )
*
*" << INTGR : 時間積分 >>
*
        CALL TINTGR             !" スペクトル時間積分
     M         ( WDTOR , WDPOR , WDDPOR, WDT   , 
     I           WTTOR , WTPOR , WTT   , 
     I           DELT  , ITRUNG, 
     C           FLAPLA, 
     C           ARAD  , DRAD             )
*
        CALL GENGD              !" 新格子点値の生成
     O         ( GAU   , GAV   , GAW   , GAT   , 
     O           GATOR , GAPOR ,
     I           WDTOR , WDPOR , WDT   , 
     C           FLAPLA, UVFACT, 
     C           ARAD  , DRAD                    )
*
      IF ( ITRUNG . LT. IORUNG ) GOTO 8888  !" <--- Runge-Kutta Loop 終わり
*
*" << OUTPUT : 出力 >>
*
      CALL HISTIN ( GTUA  , 'UA' )  !" 出力領域へ
      CALL HISTIN ( GTVA  , 'VA' )
      CALL HISTIN ( GTWA  , 'WA' )
      CALL HISTIN ( GTUT  , 'UT' )  !" 出力領域へ
      CALL HISTIN ( GTVT  , 'VT' )
      CALL HISTIN ( GTWT  , 'WT' )
      CALL HISTIN ( GTH   , 'H'  )
*
      CALL DDIAG                       !" 診断
     I         ( GAU   , GAV   , GAW   , GAT   , 
     C           PRND  , RAYL  , 
     C           ALAT  , ARAD                     )
*
      IF ( ODIASP ) THEN
        CALL SPCTLD                    !" スペクトルでの診断
     I         ( WDTOR , 'TOR' , KMAX  , NMO   )
        CALL SPCTLD
     I         ( WDPOR , 'POR' , KMAX  , NMO   )
        CALL SPCTLD
     I         ( WDT   , 'TEMP', KMAX  , NMO   )
      ENDIF
*
      RETURN
*======================================================================
      ENTRY DSTRT                      !" 初期データの再生成
     M         ( GAU   , GAV   , GAW   , GAT   , 
     O           GATOR , GAPOR ,
     C           ALAT  , DLAT  , 
     C           ARAD  , DRAD                      )
*
*   [MODIFY]
*     REAL       GAU   ( IDIM, JDIM, 0:KDIM )  !" 西風   ｕ
*     REAL       GAV   ( IDIM, JDIM, 0:KDIM )  !" 南風   ｖ
*     REAL       GAW   ( IDIM, JDIM, 0:KDIM )  !" 鉛直風 ｗ
*     REAL       GAT   ( IDIM, JDIM, 0:KDIM )  !" 温度   Ｔ
*
*   [OUTNPUT]
*     REAL       GATOR ( IDIM, JDIM, 0:KDIM )  !" トロイダル Ψ
*     REAL       GAPOR ( IDIM, JDIM, 0:KDIM )  !" ポロイダル Φ
*
*   [INPUT]
*     REAL       ALAT  ( JDIM )              !" 緯度
*     REAL       DLAT  ( JDIM )              !" 緯度荷重
*
*   [ONCE]
      IF ( .NOT. OSETC ) THEN
         OSETC = .TRUE.
         CALL    DSETC                 !" スペクトル定数設定
     O         ( NMO   , FLAPLA, EDEL  , UVFACT, 
     C           ALAT  , DLAT                      )
      ENDIF
*
      CALL STRTV                       !" データ再生成
     M         ( GAU   , GAV   , GAW   , GAT   , 
     O           GATOR , GAPOR ,
     C           FLAPLA, UVFACT,
     C           ARAD  , DRAD                     ) 
*
      RETURN
      END
***********************************************************************
      SUBROUTINE DHSTRG                !" 力学関係出力の登録
*
*   [INTERNAL PARAM]
      INTEGER    ISTYPL
      DATA       ISTYPL / 1 /
*
*   [INTERNAL WORK]
      REAL       VMISS
*
      CALL GZDBGT ( 'MISS' , VMISS )   !" 欠損値
*
      CALL HISTRG                      !" 出力登録
     I       ( 'UA', 'x-mom.advection' ,' m/s**2   ', 'ALEV',
     I           VMISS , VMISS , 2.5   , 10.   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      CALL HISTRG                      !" 出力登録
     I       ( 'VA', 'y-mom.advection' ,' m/s**2   ', 'ALEV',
     I           VMISS , VMISS , 2.5   , 10.   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      CALL HISTRG                      !" 出力登録
     I       ( 'WA', 'z-mom.advection' ,' m/s**2   ', 'ALEV',
     I           VMISS , VMISS , 2.5   , 10.   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      CALL HISTRG                      !" 出力登録
     I       ( 'UT', 'x-adv.heat flux' ,'K m/s   ', 'ALEV',
     I           VMISS , VMISS , 2.5   , 10.   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      CALL HISTRG                      !" 出力登録
     I       ( 'VT', 'y-adv.heat flux' ,'K m/s   ', 'ALEV',
     I           VMISS , VMISS , 2.5   , 10.   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      CALL HISTRG                      !" 出力登録
     I       ( 'WT', 'z-adv.heat flux' ,'K m/s   ', 'ALEV',
     I           VMISS , VMISS , 2.5   , 10.   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      CALL HISTRG                      !" 出力登録
     I       ( 'H', 'heating rate    ' ,'K/s  ', 'ALEV',
     I           VMISS , VMISS , 2.5   , 10.   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE  DPARM                 !" パラメター設定
     O             ( RAYL   , TAU  , PRND  , 
     O               ONONLN                  )
*
*   [OUTPUT]
      REAL       RAYL                          !" レイリー数
      REAL       TAU                           !" √テイラー数
      REAL       PRND                          !" プランドル数
      LOGICAL    ONONLN                        !" 非線形力学スイッチ
*
*   [INTERNAL WORK]
      REAL       RAYLI                         !" レイリー数
      REAL       TAYLR                         !" テイラー数
      REAL       PRNDL                         !" プランドル数
      LOGICAL    NONLIN                        !" 非線形力学スイッチ
*
      NAMELIST    /NMDYN/  RAYLI, TAYLR, PRNDL, NONLIN
*
      INTEGER    IFPAR, JFPAR
*
*"         < 1. NAMELIST の読み込み >
*
      RAYLI = 0.
      TAYLR = 0.
      PRNDL = 1.0
      NONLIN = .TRUE.
*
      CALL   REWNML ( IFPAR , JFPAR )
      READ   ( IFPAR, NMDYN , END=1190 )
 1190 WRITE  ( JFPAR, NMDYN  )
*
      RAYL   = RAYLI
      TAU    = SQRT( TAYLR )
      PRND   = PRNDL
      ONONLN = NONLIN
*
      RETURN
      END
