* PACKAGE DADMN  !" 力学 メインパート
*"                  ver 1.00                   92/11/16 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
*
      LOGICAL    OSETC , OFIRST
      DATA       OSETC , OFIRST  / .FALSE., .TRUE.  /
      SAVE       OSETC , OFIRST
*
*   [INTERNAL PARAM]
      LOGICAL    ODIASP
c$$$      DATA       ODIASP         / .FALSE. /
      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       GTH   ( IDIM, JDIM, 0:KDIM )  !" 温度時間変化項  Ｈ
*"       : 格子点変化項を足し込む
*
      REAL       GTUT  ( IDIM, JDIM, 0:KDIM )  !" 温度東西移流項  ＵＴ
      REAL       GTVT  ( IDIM, JDIM, 0:KDIM )  !" 温度南北移流項  ＶＴ
      REAL       GTWT  ( 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 )
*
      integer m,l,k
*
*
*"  << 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 )
         CALL    DHSTRG                    !" 出力設定
*
         CALL GD2WD             !" データをスペクトルに
     O         ( WDTOR , WDPOR , WDT   , 
     I           GATOR , GAPOR , GAT      )
*
         CALL RESET ( WDDPOR, NMDIM*KMAX )
         CALL DLOPER
     M           ( WDDPOR, 
     I             WDPOR , 
     C             FLAPLA,
     C             ARAD  , DRAD    )
      ENDIF
*
*" << Runge-Kutta Loop >>
*
      ITRUNG = 0
*
 8888 CONTINUE                    !" <---  Runge-Kutta Loop : ここから
        ITRUNG = ITRUNG + 1 
*
*" << TNDNCY : 時間変化項 >>
*
c$$$        CALL GRDDYN             !"  非線型力学項
c$$$     O         ( GTUA  , GTVA  , GTWA  , 
c$$$     O           GTUT  , GTVT  , GTWT  , GTH   ,
c$$$     I           GAU   , GAV   , GAW   , GAT   ,
c$$$     C           UVFACT, FLAPLA, 
c$$$     C           ALAT  , DLAT  , ARAD  , DRAD    )
c$$$*
c$$$        CALL TENG2W             !" 時間変化項をスペクトルに
c$$$     O         ( WTTOR , WTPOR , WTT   , 
c$$$     I           GTUA  , GTVA  , GTWA  , 
c$$$     I           GTUT  , GTVT  , GTWT  , GTH , 
c$$$     C           FLAPLA,
c$$$     C           ARAD  , DRAD                  )
        CALL RESET( WTTOR, NMDIM*KMAX )
        CALL RESET( WTPOR, NMDIM*KMAX )
        CALL RESET( WTT  , NMDIM*KMAX )
*
        CALL LNDYN              !" 力学 線型項
     M         ( WTTOR , WTPOR , WTT   , 
     I           WDTOR , WDPOR , WDT   , WDDPOR, 
     C           RAYL  , TAU   , PRND  , 
     C           NMO   , FLAPLA, 
     C           ARAD  , DRAD                    )
*
      DO 111 K = 0, KDIM
      DO 111 M = 0, MMAX
         DO 111 L = 0, LMAX
            WRITE ( 6, * )  M,L
            if( nmo( 1,m,l) .ne. 0 )then
            WRITE ( 6, * ) 
     &      WDTOR( NMO(1,M,L), K ), WTTOR( NMO(1,M,L), K )
            WRITE ( 6, * ) 
     &      WDTOR( NMO(2,M,L), K ), WTTOR( NMO(2,M,L), K )
            endif
 111     CONTINUE 
*
*" << 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 : 出力 >>
*
c$$$      CALL HISTIN ( GSIGD , 'SIGDOT')  !" 出力領域へ
c$$$      CALL HISTIN ( GPID  , 'PIDOT' )
*
      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   , 'TMP' , 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 )   !" 欠損値
*
c$$$      CALL HISTRG                      !" 出力登録
c$$$     I       ( 'UT', 'x-convective heat flux' ,'K m/s   ', 'AMLEV',
c$$$     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
c$$$     I         '   '   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
c$$$*
c$$$      CALL HISTRG
c$$$     I       ( 'VT', 'pi-tendency         ' ,'1/s   ', 'ASFC',
c$$$     I           VMISS , VMISS , VMISS , VMISS , ISTYPL  ,
c$$$     I         '   '   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE  DPARM                 !" パラメター設定
     O             ( RAYL , TAU  , PRND )
*
*   [OUTPUT]
      REAL       RAYL                          !" レイリー数
      REAL       TAU                           !" √テイラー数
      REAL       PRND                          !" プランドル数
*
*   [INTERNAL WORK]
      REAL       RAYLI                         !" レイリー数
      REAL       TAYLR                         !" テイラー数
      REAL       PRNDL                         !" プランドル数
*
      NAMELIST    /NMDYN/  RAYLI, TAYLR, PRNDL
*
      INTEGER    IFPAR, JFPAR
*
*"         < 1. NAMELIST の読み込み >
*
      RAYLI = 0.
      TAYLR = 0.
      PRNDL = 1.0
*
      CALL   REWNML ( IFPAR , JFPAR )
      READ   ( IFPAR, NMDYN , END=1190 )
 1190 WRITE  ( JFPAR, NMDYN  )
*
      RAYL = RAYLI
      TAU  = SQRT( TAYLR )
      PRND = PRNDL
*
      RETURN
      END
