* PACKAGE ALANAL !" 線型解析メイン
*                  ver 1.2                     93/09/27 takepiro
***********************************************************************
      PROGRAM GLANAL                           !" メインルーチン
*
*"    Boussinesq 全球モデル
*"      線型解析
*
*"   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZLDIM)                      !" パラメター数, 係数行列
      INCLUDE   (ZHDIM)                      !" 文字列文字数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zldim.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       GROW                  !" 成長率
      REAL       FREQ                  !" 振動数
*
      INTEGER    ISTEP                 !" 通しステップ数
      INTEGER    ITA                   !" 通し時間(t), 標準時間単位
      LOGICAL    ORSTR                 !" 固有関数出力するか否か
      INTEGER    MODE                  !" 出力固有関数のモード番号
*                                      !" ( 1 で最大成長率 )
*"       : 時刻等
*
*   [ONCE]
*"       : 実験管理パラメータ
*
*"       : 時間ステップパラメータ
*
*"    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 )      !" Δｒ(整数)
*"       : 座標値
*
*   [WORK]
      COMMON    /COMWRK/ WORK
      REAL       WORK  ( MATDMX*MATDMX ) !" ワーク領域
*
      INTEGER    MFIX( MDIM )          !" 東西波数
*
      REAL       RAYL( RDIM )          !" レイリー数
      REAL       TAYL( TDIM )          !" テイラー数
      REAL       TAU ( TDIM )          !" √テイラー数
      REAL       PRND( PDIM )          !" プランドル数
*
      INTEGER    MSTEP , MEND          !" 東西波数の計算数
      INTEGER    RSTEP , REND          !" レイリー数の計算数
      INTEGER    TSTEP , TEND          !" テイラー数の計算数
      INTEGER    PSTEP , PEND          !" プランドル数の計算数
*
      INTEGER    IUNIT                 !" 標準出力装置番号
*
#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 GLPGET( 'MSGUNIT', IUNIT )
      WRITE ( IUNIT, * )' LINEAR ANALYSIS GLANAL'
     &                   //NAME
     &                   //' ver.1.2, 93/09/27'
*
      CALL CLCSTR ( 'SETUP' )
      CALL YPREP                       !" システム前処理 
*
      CALL SETLPR                      !" パラメター範囲設定
     O       ( MODE   ,
     O         MFIX   , MEND   ,
     O         RAYL   , TAYL   , TAU    , PRND   ,
     O         REND   , TEND   , PEND              )
*
      CALL SETCOR                      !" 座標値
     O         ( ALON  , DLON  ,
     O           ALAT  , DLAT  ,
     O           ARAD  , DRAD    )
*
c$$$      CALL AHSTRG                      !" 標準時間平均出力の登録
*
      CALL CLCEND ( 'SETUP' )
*
*  << LOOP : <************* メイン・ループ >>
*
      DO 9000 MSTEP = 1, MEND
         DO 8000 PSTEP = 1, PEND
            DO 8000 TSTEP = 1, TEND
               DO 8000 RSTEP = 1, REND
*
*     << START : ステップはじめの処理 >>
*
                  CALL CLCSTR( 'SETLTM' )
*
                  CALL SETLTM
     O                   (  ITA   , ISTEP , ORSTR ,
     I                      MSTEP , RSTEP , TSTEP , PSTEP  )
*
                  CALL SETTIM                       !" 波数を時刻として登録
     I                   ( MFIX(MSTEP)   ,  ISTEP  )
*
                  CALL HSPARM
     I               ( 'WAVE NUMBER'     , REAL(MFIX(MSTEP)) )
                  CALL HSPARM
     I               ( 'RAYLEIGH NUMBER' , RAYL(RSTEP)    )
                  CALL HSPARM
     I               ( 'TAYLOR NUMBER'   , TAU(TSTEP)**2  )
                  CALL HSPARM
     I               ( 'PRANDTL NUMBER'  , PRND(PSTEP)    )
*
                  WRITE ( IUNIT, * ) ' m  : ', MFIX(MSTEP) 
                  WRITE ( IUNIT, * ) ' Ra : ', RAYL(RSTEP) 
                  WRITE ( IUNIT, * ) ' Ta : ', TAU(TSTEP)**2
                  WRITE ( IUNIT, * ) ' Pr : ', PRND(PSTEP)
*
                  CALL CLCEND ( 'SETLTM' )
*

*"    << INTEGR : 時間積分 >>
* 
c$$$                   CALL CLCSTR ( 'DLANAL' )
*
                   CALL LANAL
     O                   ( GAU   , GAV   , GAW   , GAT   ,
     O                     GATOR , GAPOR , 
     O                     GROW  , FREQ  , 
     I                     MFIX(MSTEP), 
     I                     RAYL(RSTEP), TAU(TSTEP), PRND(PSTEP), 
     I                     ORSTR , MODE  , 
     C                     ALON  , DLON  , ALAT  , DLAT  ,
     C                     ARAD  , DRAD                    )
*
c$$$                   CALL CLCEND ( 'DLANAL' )
*
*"    << OUTPUT : データ出力 >>
*
                   CALL CLCSTR( 'OUTPUT' )
*
                   CALL LHSEIN
     I                  (  GROW  , FREQ   , 
     I                     MFIX   , 
     I                     RSTEP  , TSTEP  , PSTEP  ) 
                   CALL HSPARM
     I                  ( 'GROWTH RATE', GROW )
                   CALL HSPARM
     I                  ( 'FREQENCY'   , FREQ  )
*
                   IF ( ORSTR ) THEN
                      WRITE ( IUNIT, * ) ' GROWTH RATE : ', GROW
                      WRITE ( IUNIT, * ) ' FREQENCY    : ', FREQ
*
                      CALL ANORML
     I                     ( GAU   , GAV   , GAW   , GAT   , 
     I                       GATOR , GAPOR , 
     I                       'GAT'                           )
*
                      CALL WRRSTR      !" リスタート・ファイルの書き込み
     I                     ( GAU   , GAV   , GAW   , GAT   , 
     I                       GATOR , GAPOR ,
     I                       MFIX(MSTEP)   , ISTEP              )
                   ENDIF
*
                   CALL CLCEND( 'OUTPUT' )
*
*" << LEND : メイン・ループ終わり *************> >>
*
 8000    CONTINUE 
*
         CALL LHSEOU
     I          (  MFIX(MSTEP)  , 
     C             REND   , TEND   , PEND   ) 
*
 9000 CONTINUE 
*
      CALL CLCOUT                      !" CPU時間出力
      CALL YFINE                       !" システム後処理
*
      STOP
      END
***********************************************************************
      SUBROUTINE AHSTRG                !" 標準時間平均出力の登録
*
*   [INTERN PARAM]
      INTEGER    ISTYPL
      DATA       ISTYPL / 1 /
      REAL       VMISS
*
      CALL GZDBGT ( 'MISS', VMISS )   !" 欠損値
*
      CALL     HISTRG                  !" 出力の登録
     I       ( 'U     ', 'u-velocity          ' ,'m/s   ', 'ALEV',
     I           VMISS , VMISS , 0.1   , 0.5   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      CALL     HISTRG
     I       ( 'V     ', 'v-velocity          ' ,'m/s   ', 'ALEV',
     I           VMISS , VMISS , 0.1   , 0.5   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      CALL     HISTRG
     I       ( 'W     ', 'w-velocity          ' ,'m/s   ', 'ALEV',
     I           VMISS , VMISS , 0.1   , 0.5   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      CALL     HISTRG
     I       ( 'T     ', 'temperature         ' ,'K     ', 'ALEV',
     I           VMISS , VMISS , 0.1   , 0.5   , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , ' '     ,'(F12.3)' )
*
      CALL     HISTRG
     I       ( 'TOR   ', 'toroidal potential  ' ,'1/s   ', 'ALEV',
     I           VMISS , VMISS , 0.1   , 0.5   , ISTYPL  ,
     I         '   '   , '   ' , 0     , 0     , ' '     ,'(1PE12.3)' )
*
      CALL     HISTRG
     I       ( 'POR   ', 'poroidal potential  ' ,'1/s   ', 'ALEV',
     I           VMISS , VMISS , 0.1   , 0.5   , ISTYPL  ,
     I         '   '   , '   ' , 0     , 0     , ' '     ,'(1PE12.3)' )
*
      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
**********************************************************************
      SUBROUTINE ANORML
     M            ( GAU   , GAV   , GAW   , GAT   , 
     M              GATOR , GAPOR , 
     I              HITEM                            )
*
*"   [PARAM]
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#endif
*
*   [MODIFY]
      REAL       GAU   ( IDIM, JDIM, KMAX )  !" 西風   ｕ
      REAL       GAV   ( IDIM, JDIM, KMAX )  !" 南風   ｖ
      REAL       GAW   ( IDIM, JDIM, KMAX )  !" 鉛直風 ｗ
      REAL       GAT   ( IDIM, JDIM, KMAX )  !" 温度   Ｔ
      REAL       GATOR ( IDIM, JDIM, KMAX )  !" トロイダル Ψ
      REAL       GAPOR ( IDIM, JDIM, KMAX )  !" ポロイダル Φ
*
*   [INPUT]
      CHARACTER  HITEM * 5                    !" 規格化変数名
*
*   [INTERNAL WORK]
      REAL       ANORM  , ABSMAX
      INTEGER    I, J, K
*
      IF( HITEM(1:3) .EQ. 'GATOR' )THEN
         ANORM = ABSMAX
     I         ( GATOR, IDIM, JDIM, KMAX, IMAX, JMAX )
      ELSE IF( HITEM(1:3) .EQ. 'GAPOR' )THEN
         ANORM = ABSMAX
     I         ( GAPOR, IDIM, JDIM, KMAX, IMAX, JMAX )
      ELSE IF( HITEM(1:3) .EQ. 'GAU' )THEN
         ANORM = ABSMAX
     I         ( GAU, IDIM, JDIM, KMAX, IMAX, JMAX )
      ELSE IF( HITEM(1:3) .EQ. 'GAV' )THEN
         ANORM = ABSMAX
     I         ( GAV, IDIM, JDIM, KMAX, IMAX, JMAX )
      ELSE IF( HITEM(1:3) .EQ. 'GAW' )THEN
         ANORM = ABSMAX
     I         ( GAW, IDIM, JDIM, KMAX, IMAX, JMAX )
      ELSE IF( HITEM(1:3) .EQ. 'GAT' )THEN
         ANORM = ABSMAX
     I         ( GAT, IDIM, JDIM, KMAX, IMAX, JMAX )
      ELSE
         RETURN
      ENDIF
*
      DO 8100 K=1,KMAX
         DO 8100 J=1,JMAX
            DO 8100 I=1,IMAX
               GAU(I,J,K) = GAU(I,J,K) / ANORM 
               GAV(I,J,K) = GAV(I,J,K) / ANORM 
               GAW(I,J,K) = GAW(I,J,K) / ANORM 
               GAT(I,J,K) = GAT(I,J,K) / ANORM 
               GATOR(I,J,K) = GATOR(I,J,K) / ANORM 
               GAPOR(I,J,K) = GAPOR(I,J,K) / ANORM 
 8100 CONTINUE
      RETURN
      END
**********************************************************************
      FUNCTION ABSMAX      !"  配列の絶対値の最大値を返す
     I         ( DATA, IDIM, JDIM, KMAX, IMAX, JMAX )
*
*   [PARAM] 
      INTEGER    IDIM
      INTEGER    JDIM
      INTEGER    KMAX
*
*   [INPUT] 
      REAL       DATA( IDIM, JDIM, KMAX )
      INTEGER    IMAX
      INTEGER    JMAX
*
*   [OUTPUT] 
      REAL       ABSMAX
*
*   [INTERNAL WORK] 
      INTEGER    I, J, K
*
      ABSMAX = -1.E36
*
      DO 1100 K=1,KMAX
        DO 1100 J=1,JMAX
          DO 1100 I=1,IMAX
            IF ( ABS(DATA( I,J,K )) .GT. ABSMAX ) THEN
              ABSMAX = DATA( I,J,K )
            ENDIF
 1100 CONTINUE
*
      RETURN
      END
