*" PACKAGE AEDYFL !" 渦フラックス計算
*"                  ver 1.00                   93/06/25 takepiro
***********************************************************************
      PROGRAM AEDYFL                        !" メインルーチン
*
*"    Boussinesq 全球モデル
*"      非線形項の計算, 出力
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZHDIM)                      !" 文字列文字数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zhdim.F"                    !" 文字列文字数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [VAR]
      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)
*
      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 )  !" 温度時間変化項  Ｈ
*"       : <DYNMCS> 格子点変化項(D)
*
      REAL       WTUA  ( NMDIM , KMAX     )  !" 
      REAL       WDDUA ( NMDIM , KMAX     )  !" 
      REAL       CDATA ( NMDIM , KMAX     )  !" 
*
      INTEGER    ISTEP                 !" 通しステップ数
      INTEGER    ITA                   !" 通し時間(t), 標準時間単位
c$$$      INTEGER    IDATEA( 6 )           !" 時刻年月日時分秒(ダミー)
*"       : 時刻等
*
*   [ONCE]
      INTEGER    ITSTRT                !" 計算開始時刻(標準時間単位)
      INTEGER    ITEND                 !" 計算終了時刻(標準時間単位)
      INTEGER    IORSTR                !" 出力間隔：再出発(標準時間単位)
*"       : 実験管理パラメータ
*
      INTEGER    IDELT                 !" 標準時間刻み(SEC)
*"       : 時間ステップパラメータ
*
*"    COMMON /COMCON/ ( include ZCCOM )
*     REAL       ERI                          !" 内径
*     REAL       ERO                          !" 外径
*     REAL       ETA                          !" 内径/外径
*"    COMMON /COMCON/ end
*"       : 物理パラメータ
*
      REAL       ALAT  ( JDIM )        !" 緯度
      REAL       DLAT  ( JDIM )        !" 緯度荷重
      REAL       ALON  ( IDIM )        !" 経度
      REAL       DLON  ( IDIM )        !" 経度荷重
      REAL       ARAD  ( 0:KDIM )        !" ｒレベル(整数)
      REAL       DRAD  ( 0:KDIM )        !" Δｒ(整数)
*"       : 座標値
*
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       EDEL  ( NMDIM  )            !" ζ，Ｄ→Ｕ，Ｖ
      REAL       UVFACT ( IDIM, JDIM )       !" u→U のファクター
*
*   [WORK]
      COMMON    /COMWRK/ WORK
      REAL       WORK  ( NWORK )       !" ワーク領域
*
c$$$      DATA      IDATEA / 1992, 12, 18, 0, 0, 0/
*
      INTEGER  K, L, NM
*
#ifdef DEBUG
      external common_handler
      integer i, ieee_handler
      i = ieee_handler("set","common",common_handler )
      if( i.ne. 0 ) print *, "Could not establish fp signal handler"
#endif
*
*" << SETPUP : 初期設定 >>
*
      CALL CLCSTR ( 'SETUP' )
      CALL YPREP                       !" システム前処理 
*
      CALL SETPAR                      !" 実験パラメータ
     O         ( ITSTRT , ITEND , IORSTR ,
     O           IDELT                      )
*
      CALL SETCOR                      !" 座標値
     O         ( ALON  , DLON  ,
     O           ALAT  , DLAT  ,
     O           ARAD  , DRAD   )
*
      CALL RDSTRT                      !" 初期値の読み込みと再生成
     O         ( GAU   , GAV   , GAW   , GAT   , 
     O           GATOR , GAPOR ,
     O           ITA   , ISTEP ,
     C           ITSTRT, ALAT  , DLAT  , ARAD  , DRAD     )
*
      CALL ADMXMN 
     I         ( GAU   , GAV   , GAW   , GAT   , 
     I           GATOR , GAPOR , ITA   , 'start GA'      )
*
      CALL DHSTRG                      !" 出力設定

      CALL SETTIM                      !" 時刻を記憶 
     I         ( ITA-IDELT, ISTEP-1  )
*
c$$$      CALL HISTOU
*
      CALL SETTIM                      !" 時刻を記憶 
     I         ( ITA   , ISTEP  )
*
*
*"  << 非線形項の計算 >>
*
      CALL DSETC                       !" スペクトル定数設定
     O         ( NMO   , FLAPLA, EDEL  , UVFACT, 
     C           ALAT  , DLAT                      )
*
      CALL GRDDYN                      !"  非線型力学項
     O      ( GTUA  , GTVA  , GTWA  , 
     O        GTUT  , GTVT  , GTWT  , GTH   ,
     I        GAU   , GAV   , GAW   , GAT   ,
     C        UVFACT, FLAPLA, 
     C        ALAT  , DLAT  , ARAD  , DRAD    )
*
      CALL G2W
     O         ( WTUA ,
     I           GTUA  ,
     I           '  '  , 'POS ',  KMAX  )
*
      DO 3000 K = 1, KMAX
         DO 3000 NM = 1, NMDIM
            WDDUA ( NM, K ) = WTUA ( NM, K ) / ARAD( K )
 3000 CONTINUE 
*
      CALL W2C
     M         ( CDATA ,
     I           WTUA ,
     I           ' '   , 'POS ', KMAX )
*
      CALL C2W
     M         ( WDDUA ,
     I           CDATA ,
     I           'DR ' , 'ADD ', KMAX )
*
      DO 3000 K = 1, KMAX
         DO 3000 M = 0, MMAX 
            DO 3000 L = 0, LMAX
               DO 3000 I = 1, 2
                  NM = NMO( I, M, L )
                  IF ( NM .NE. 0 )THEN
                     WRITE( 6, * ) M,L, K, WDDUA( NM,K )
                     ENDIF
                  ENDIF
 3000 CONTINUE 
*
*"  << 非線形項の出力 >>
*
c$$$      CALL WREDFL            !" 渦フラックスのファイル出力
c$$$     I      ( GTUA  , GTVA  , GTWA  , 
c$$$     I        GTUT  , GTVT  , GTWT  , 
c$$$     I        GTH                      )
*
c$$$      CALL HISTIN ( GTUA  , 'UA' )  !" 出力領域へ
c$$$      CALL HISTIN ( GTVA  , 'VA' )
c$$$      CALL HISTIN ( GTWA  , 'WA' )
c$$$      CALL HISTIN ( GTUT  , 'UT' )  !" 出力領域へ
c$$$      CALL HISTIN ( GTVT  , 'VT' )
c$$$      CALL HISTIN ( GTWT  , 'WT' )
c$$$      CALL HISTIN ( GTH   , 'H'  )
c$$$*
c$$$      CALL HISTOU
*
      CALL YFINE                       !" システム後処理
*
      STOP
      END
***********************************************************************
      SUBROUTINE WREDFL            !" 渦フラックスのファイル出力
     I         ( GTUA  , GTVA  , GTWA  , 
     I           GTUT  , GTVT  , GTWT  , 
     I           GTH                      )

*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZHDIM)                      !" 文字列文字数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zhdim.F"                    !" 文字列文字数
#endif
*
*   [INPUT]
      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 )  !" 温度時間変化項  Ｈ
*"       : 格子点データ(t) <DYNMCS>  格子点データ(t+Δt)
*
      INTEGER    ISTEP                 !" 通しステップ数
      INTEGER    ITA                   !" 通し時間(t)
      INTEGER    ITB                   !" 通し時間(t-Δt)
*"       : 時刻等
*
*    [INTERNAL PARAM]
      INTEGER    IDATEA( 6 )           !" 時刻：年月日時分秒
      INTEGER    IDATEB( 6 )           !" 時刻：年月日時分秒
      CHARACTER  FILE   *(NFILN)       !" 出力ファイル：再出発
      CHARACTER  DFMT  *(NCC)          !" 出力フォーマット
      CHARACTER  AVRTYP*(NCC)          !" 平均種類
      DATA       FILE   / '  ' /
      DATA       DFMT   / 'UR8'/
      DATA       AVRTYP / ' '  /
      INTEGER    JFILED
      DATA       JFILED / 41 /
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      LOGICAL    OADVNC
      REAL       DELT
*
      INTEGER    IOMODE
      INTEGER    NOEND
      INTEGER    ITDUR
      DATA       IOMODE / 0  /
      DATA       NOEND  / 0  /
      DATA       ITDUR  / 1  /
*
      SAVE
*
      NAMELIST  /NMFEDY/ FILE, DFMT, AVRTYP
*    
*    [INTERNAL WORK]
      INTEGER    IFPAR, JFPAR, JFLRST, IERR
*
      IF ( OFIRST ) THEN
         OFIRST = .FALSE.
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMFEDY, END=1090 )
 1090    WRITE  ( JFPAR, NMFEDY )         
*
         CALL IFLOPN
     O         ( JFLRST, IERR  ,
     I           FILE  , JFILED, 'WRITE', 'UNFORMATTED' )
         CALL INQTIM
     O         ( ITA   , ITB   , IDATEA, IDATEB,
     O           ISTEP , DELT  , OADVNC          )
      ENDIF
*
*" << GAWRIT : ファイル出力 (GA) >>
*
      CALL EDYWRT
     I        ( GTUA  ,
     I         'GTUA' , 'x-mom.adv. @ t       ','m/s**2' , AVRTYP , 
     I          ITA   , IDATEA, ISTEP , ITDUR ,
     I          JFLRST, IOMODE, NOEND , 'ALEV', DFMT             )
*
      CALL EDYWRT
     I        ( GTVA  ,
     I         'GTVA' , 'y-mom.adv. @ t       ','m/s**2' ,AVRTYP , 
     I          ITA   , IDATEA, ISTEP , ITDUR ,
     I          JFLRST, IOMODE, NOEND , 'ALEV', DFMT            )
*
      CALL EDYWRT
     I        ( GTWA  ,
     I         'GTWA' , 'z-mom.adv. @ t       ','m/s**2' , AVRTYP , 
     I          ITA   , IDATEA, ISTEP , ITDUR ,
     I          JFLRST, IOMODE, NOEND , 'ALEV', DFMT             )
*
      CALL EDYWRT
     I        ( GTUT  ,
     I         'GTUT' , 'x-adv.heat flux @ t  ','K m/s' , AVRTYP , 
     I          ITA   , IDATEA, ISTEP , ITDUR ,
     I          JFLRST, IOMODE, NOEND , 'ALEV', DFMT            )
*
      CALL EDYWRT
     I        ( GTVT  ,
     I         'GTVT' , 'y-adv.heat flux @ t  ','K m/s' , AVRTYP , 
     I          ITA   , IDATEA, ISTEP , ITDUR ,
     I          JFLRST, IOMODE, NOEND , 'ALEV', DFMT           )
*
      CALL EDYWRT
     I        ( GTWT  ,
     I         'GTWT' , 'z-adv.heat flux @ t  ','K m/s' , AVRTYP , 
     I          ITA   , IDATEA, ISTEP , ITDUR ,
     I          JFLRST, IOMODE, NOEND , 'ALEV', DFMT            )
*
      CALL EDYWRT
     I        ( GTH   ,
     I         'GAH'  , 'heating rate @ t    ','K/s' , AVRTYP , 
     I          ITA   , IDATEA, ISTEP , ITDUR ,
     I          JFLRST, IOMODE, NOEND , 'ALEV', DFMT         )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE EDYWRT      !" GTOOL3出力
     I         ( GDATA ,
     I           HITEM , HTITL , HUNIT , HAVTYP,
     I           IT    , IDATE , ISTEP , ITDUR ,
     I           JFILE , IOMODE, NOEND , HCORD , HDFMT  )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZHDIM)                      !" 文字数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zhdim.F"                    !" 文字数
#endif
*
*    [INPUT]     
      REAL       GDATA ( IDIM, JDIM, * )     !" データ
      CHARACTER  HITEM  *(*)                 !" 識別名称
      CHARACTER  HTITL  *(*)                 !" タイトル
      CHARACTER  HUNIT  *(*)                 !" 単位
      CHARACTER  HAVTYP *(*)                 !" 平均種類
      INTEGER    IT                          !" 通し時間
      INTEGER    IDATE ( * )                 !" 時刻：年月日時分秒
      INTEGER    ISTEP                       !" 通しステップ数
      INTEGER    ITDUR                       !" 代表時間
      INTEGER    JFILE                       !" 出力ファイル番号
      INTEGER    IOMODE                      !" 出力モード
      INTEGER    NOEND                       !" ENDFILE間隔
      CHARACTER  HCORD  *(*)                 !" 座標系
      CHARACTER  HDFMT        *(NCC)         !" データフォーマット
*
*    [INTERNAL WORK]
      REAL       GDATAW ( IDIM, JDIM, KMAX )
      CHARACTER  HHEAD ( NDC )*(NCC)         !" ヘッダー
      REAL       VMISS
      INTEGER    ISIZE
*
      CALL GZDBGT( 'MISS' , VMISS )
*
      CALL RESET ( GDATAW, IDIM*KDIM*KMAX )
*
      CALL HSTADD
     M         ( GDATAW,
     I           GDATA , 1      ,
     I           HCORD , HAVTYP  )
*
      CALL EDYHED                    !" ヘッダーの設定
     O         ( HHEAD , ISIZE ,
     I           HCORD , HAVTYP,
     I           HITEM , HTITL , HUNIT ,
     I           HDFMT , VMISS ,
     I           VMISS , VMISS , VMISS , VMISS , 1      )
*
      CALL HHPARM                    !" パラメターのヘッダー登録
     M      ( HHEAD )
*
      CALL GZDBWT                    !" 精度変換書き込み
     O            ( HHEAD , GDATAW,
     O              IT    , IDATE , ISTEP , ITDUR ,
     I              JFILE , IOMODE, NOEND , HDFMT   )
*
      RETURN
      END
***********************************************************************
      SUBROUTINE EDYHED       !" ヘッダーのセット
     O         ( HHEAD , ISIZE ,
     I           HCORD , HAVTYP,
     I           HITEM , HTITL , HUNIT ,
     I           HDFMT , VMISS ,
     I           DMIN  , DMAX  , DIVS  , DIVL  , ISTYP  )
*
*   [OUTPUT] 
      CHARACTER  HHEAD ( * )*(*)             !" ヘッダー
      INTEGER    ISIZE                       !" サイズ
*
*   [INPUT] 
      CHARACTER  HCORD  *(*)                 !" 座標系
      CHARACTER  HAVTYP *(*)                 !" 平均種類
      CHARACTER  HITEM  *(*)                 !" 識別名称
      CHARACTER  HTITL  *(*)                 !" タイトル
      CHARACTER  HUNIT  *(*)                 !" 単位
      CHARACTER  HDFMT  *(*)                 !" データフォーマット
      REAL       VMISS                       !" 欠損値の値
      REAL       DMIN                        !" レンジ(最小)
      REAL       DMAX                        !" レンジ(最大)
      REAL       DIVS                        !" 間隔(小)
      REAL       DIVL                        !" 間隔(大)
      INTEGER    ISTYP                       !" スケーリングタイプ
*
*   [INTERNAL PARM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZHDIM)                      !" 文字列文字数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zhdim.F"                    !" 文字列文字数
#endif
*   [INTERNAL WORK] 
      CHARACTER  HAITMX *(NCC)               !" 第1軸名称
      CHARACTER  HAITMY *(NCC)               !" 第2軸名称
      CHARACTER  HAITMZ *(NCC)               !" 第3軸名称
      CHARACTER  HUTIM  *(NCC)               !" 時間単位
      CHARACTER  HEDIT  *(NCC)               !" 時刻単位
      CHARACTER  HETTL  *(NCC)               !" 時刻単位
*
      INTEGER    IXDIM, IYDIM, IZDIM, KLEVS
      INTEGER    NAVR
*
*   [EXTERNAL FUNC]
      INTEGER    LENC
*
*"         < 1. 軸名称，大きさ >
*
      CALL HQCORN
     I     ( HCORD ,
     O       HAITMX, HAITMY, HAITMZ,
     O       KLEVS                   )
*
      IF      ( HAVTYP .EQ. 'X' ) THEN
         HAITMX = ' '
         HAITMY = HAITMY
         HAITMZ = HAITMZ
         IXDIM  = 1
         IYDIM  = JMAX
         IZDIM  = KLEVS
      ELSE IF ( HAVTYP .EQ. 'Z' ) THEN
         HAITMZ = ' '
         IXDIM  = IMAX
         IYDIM  = JMAX
         IZDIM  = 1
      ELSE IF ( HAVTYP .EQ. 'XZ' ) THEN
         HAITMX = HAITMY
         HAITMY = ' '
         HAITMZ = ' '
         IXDIM  = JMAX
         IYDIM  = 1
         IZDIM  = 1
      ELSE IF ( HAVTYP .EQ. 'XY' ) THEN
         HAITMX = HAITMZ
         HAITMY = ' '
         HAITMZ = ' '
         IXDIM  = KLEVS
         IYDIM  = 1
         IZDIM  = 1
      ELSE IF ( HAVTYP .EQ. 'XYZ' ) THEN
         HAITMX = ' '
         HAITMY = ' '
         HAITMZ = ' '
         IXDIM  = 1
         IYDIM  = 1
         IZDIM  = 1
      ELSE
         IXDIM  = IMAX
         IYDIM  = JMAX
         IZDIM  = KLEVS
      ENDIF
      ISIZE = IXDIM * IYDIM * IZDIM
*
*"         < 2. ヘッダーのパック >
*
      CALL AQUTIM ( HUTIM )
*
      CALL GZDBHP                    !" 識別記述子のパック
     O         ( HHEAD ,
     I           ' '   , HITEM , 1     , 1     ,
     I           HTITL , HUNIT ,
     I           0     , ' '   , HUTIM , 1     ,
     I           HAITMX, 1     , IXDIM ,
     I           HAITMY, 1     , IYDIM ,
     I           HAITMZ, 1     , IZDIM ,
     I           HDFMT , VMISS ,
     I           DMIN  , DMAX  , DIVS  , DIVL  ,
     I           ISTYP ,
     I           ' '   , ' '                             )
*
      IF ( HAVTYP .NE. ' ' ) THEN
         NAVR  = LENC( HAVTYP )
         HEDIT = HAVTYP(1:NAVR)//'-AVG'
         HETTL = HAVTYP(1:NAVR)//'-average'
         CALL GHEADD ( HHEAD, HEDIT, HETTL )
      ENDIF
*
      RETURN
      END
***********************************************************************
      SUBROUTINE ADMXMN        !" debug monitor
     I         ( GDU   , GDV   , GDW   , GDT   , 
     I           GDTOR , GDPOR , IT    , HLABEL          )
*
*   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#endif
*
*   [INPUT] 
      REAL       GDU   ( IDIM, JDIM, KMAX )  !" 西風   ｕ
      REAL       GDV   ( IDIM, JDIM, KMAX )  !" 南風   ｖ
      REAL       GDW   ( IDIM, JDIM, KMAX )  !" 鉛直風 ｗ
      REAL       GDT   ( IDIM, JDIM, KMAX )  !" 温度   Ｔ
      REAL       GDTOR ( IDIM, JDIM, KMAX )  !" トロイダル Ψ
      REAL       GDPOR ( IDIM, JDIM, KMAX )  !" ポロイダル Φ

*
      INTEGER    IT
      CHARACTER  HLABEL *(*)
*
      WRITE ( 6,* ) HLABEL, ' IT= ', IT
      CALL MAXMIN( GDU,   IDIM, JDIM, KMAX, IDIM, JDIM, 'U'   )
      CALL MAXMIN( GDV,   IDIM, JDIM, KMAX, IDIM, JDIM, 'V'   )
      CALL MAXMIN( GDW,   IDIM, JDIM, KMAX, IDIM, JDIM, 'W'   )
      CALL MAXMIN( GDT,   IDIM, JDIM, KMAX, IDIM, JDIM, 'T'   )
      CALL MAXMIN( GDTOR, IDIM, JDIM, KMAX, IDIM, JDIM, 'TOR' )
      CALL MAXMIN( GDPOR, IDIM, JDIM, KMAX, IDIM, JDIM, 'POR' )
*
      RETURN
      END
