* PACKAGE ALCRT !" 線型解析メイン(臨界値)
*                  ver 1.2                   93/09/27 takepiro
***********************************************************************
      PROGRAM GLCRT                          !" メインルーチン
*
*"    Boussinesq 全球モデル
*"      線型解析(臨界値)  : Bracket method ( はさみうち法? )
*
*"   [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       FREQ                        !" 振動数
*
      INTEGER    ISTEP                 !" 通しステップ数
      INTEGER    ITA                   !" 通し時間(t), 標準時間単位
      LOGICAL    ORSTR                 !" 再出発出力するか否か
*"       : 時刻等
*
*   [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          !" プランドル数の計算数
*
      REAL       RAYLM                 !" レイリー数
      REAL       GRMXM                 !" 最大成長率
      INTEGER    IT    , ITMAX         !" 反復回数, 最大反復回数
      REAL       RAYERR                !" レイリー数の相対許容誤差
      REAL       RAYLL , RAYLU         !" レイリー数下限, 上限
      REAL       GRMXL , GRMXU         !" 最大成長率
*
      LOGICAL    OCONV                 !" 収束判定スイッチ
      INTEGER    MODE                  !" 出力固有関数のモード番号
*                                      !" ( 1 で最大成長率 )
      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 GLCRT'
     &                   //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    )
*
      CALL CLCEND ( 'SETUP' )
*
*  << LOOP : <************* メイン・ループ >>
*
      DO 9000 MSTEP = 1, MEND
         DO 8000 PSTEP = 1, PEND
            DO 8000 TSTEP = 1, TEND
*
*     << START : ステップはじめの処理 >>
*
               CALL CLCSTR( 'SETLTM' )
*
               CALL SETCPR                 !" 臨界値計算用パラメター設定
     O                ( ITMAX  , RAYERR , 
     O                  RAYLL  , RAYLU    )
*
               CALL SETLTM
     O                   (  ITA   , ISTEP , ORSTR ,
     I                      MSTEP , 1     , TSTEP , PSTEP  )
*
               CALL SETTIM                        !" 波数を時刻として登録
     I                   ( MFIX(MSTEP)   ,  ISTEP  )
*
               CALL HSPARM
     I               ( 'WAVE NUMBER'    ,  REAL(MFIX(MSTEP)) )
               CALL HSPARM
     I               ( 'TAYLOR NUMBER'  ,  TAU(TSTEP)**2  )
               CALL HSPARM
     I               ( 'PRANDTL NUMBER' ,  PRND(PSTEP)    )
*
               WRITE ( IUNIT, * ) ' m  : ', MFIX(MSTEP) 
               WRITE ( IUNIT, * ) ' Ta : ', TAU(TSTEP)**2
               WRITE ( IUNIT, * ) ' Pr : ', PRND(PSTEP)
*
               CALL CLCEND ( 'SETLTM' )
*
               CALL LANAL
     O                ( GAU   , GAV   , GAW   , GAT   ,
     O                  GATOR , GAPOR , 
     O                  GRMXL , FREQ , 
     I                  MFIX(MSTEP), 
     I                  RAYLL , TAU(TSTEP), PRND(PSTEP), 
     I                  ORSTR , MODE ,
     C                  ALON  , DLON  , ALAT  , DLAT  ,
     C                  ARAD  , DRAD                    )
*
               IF ( GRMXL .GT. 0 )THEN
		  WRITE( 6, * ) GRMXL
                  CALL MSGDMP
     &                   ('E','ALCRT',
     &                    'LOWER RAYLEIGH NUMBER SHOULD BE SMALLER')
               ENDIF
*
               CALL LANAL
     O                ( GAU   , GAV   , GAW   , GAT   ,
     O                  GATOR , GAPOR , 
     O                  GRMXU , FREQ  , 
     I                  MFIX(MSTEP), 
     I                  RAYLU , TAU(TSTEP), PRND(PSTEP), 
     I                  ORSTR , MODE  , 
     C                  ALON  , DLON  , ALAT  , DLAT  ,
     C                  ARAD  , DRAD                    )
*
               IF ( GRMXU .LT. 0 )THEN
                  CALL MSGDMP
     &                   ('E','ALCRT',
     &                    'UPPER RAYLEIGH NUMBER SHOULD BE LARGER')
               ENDIF
*
               RAYLM = RAYLU
               GRMXM = GRMXU
*
*"    << 分割法 loop START ------------------------------>>
* 
               DO 1000 IT = 1, ITMAX
*
c$$$                  CALL UBISEC          !" 2 分法による解の推定
c$$$                  CALL URFLSP          !" False Point 法による解の推定
                  CALL UBRENT          !" Brent 法による解の推定
     M                 ( RAYLL , RAYLU  , RAYLM  , 
     M                   GRMXL , GRMXU  , GRMXM  , 
     O                   OCONV  ,
     I                   RAYLM*RAYERR              )
*
                  IF ( OCONV ) GOTO 1100                !" iteration 終了
*
                  CALL LANAL
     O                   ( GAU   , GAV   , GAW   , GAT   ,
     O                     GATOR , GAPOR , 
     O                     GRMXM , FREQ  , 
     I                     MFIX(MSTEP), 
     I                     RAYLM , TAU(TSTEP), PRND(PSTEP), 
     I                     ORSTR , MODE  , 
     C                     ALON  , DLON  , ALAT  , DLAT  ,
     C                     ARAD  , DRAD                    )
*
 1000          CONTINUE 
*
               CALL MSGDMP                              !" 収束せず
     &                ('W','ALCRT',
     &                 'RAYLEIGH NUMBER NOT CONVERGE' )
*
*"    <<-------------------------------- 分割法 loop END >>
*
 1100          CONTINUE 
*
               CALL HSPARM
     I               ( 'RAYLEIGH NUMBER',  RAYLM  )
               CALL HSPARM
     I               ( 'GROWTH RATE',      GRMXM  )
               CALL HSPARM
     I               ( 'FREQENCY',         FREQ   )
*
*
*"    << OUTPUT : データ出力 >>
*
               CALL CLCSTR( 'OUTPUT' )
*
               CALL LHSEIN
     I                (  GRMXM  , FREQ   , 
     I                   MFIX   , 
     I                   1      , TSTEP  , PSTEP  ) 
               CALL LHSRIN
     I                (  RAYLM  , 
     I                   MFIX   , 
     I                   1      , TSTEP  , PSTEP  ) 
*
                  WRITE ( IUNIT, * ) ' RAYLEIGH NUMBER : ', RAYLM
                  WRITE ( IUNIT, * ) ' FREQENCY        : ', FREQ
                  WRITE ( IUNIT, * ) ' ITERATION       : ', IT
                  WRITE ( IUNIT, * ) ' GROWTH RATE     : ', GRMXM
                  WRITE ( IUNIT, * ) ' R. NUMBER ERROR : ', 
     &             MAX( ABS(RAYLU-RAYLM), ABS(RAYLM-RAYLM) ) / RAYLM
                  WRITE ( IUNIT, * ) 
*
               IF ( ORSTR ) THEN  !" 温度の最大値で規格化
                  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 LHSREO
     I          (  MFIX(MSTEP)  , 
     C             1    , TEND   , PEND   ) 
*
 9000 CONTINUE 
*
      CALL CLCOUT                      !" CPU時間出力
      CALL YFINE                       !" システム後処理
*
      STOP
      END
***********************************************************************
      SUBROUTINE SETCPR                !" 臨界値計算用パラメター設定
     O             ( IMAX , ERR  , 
     O               RL   , RU       )
*
      INTEGER    IMAX                  !" 最大反復回数
      REAL       ERR                   !" レイリー数の相対許容誤差
      REAL       RL    , RU            !" レイリー数下限, 上限
*
      INTEGER    ITMAX                 !" 最大反復回数
      REAL       RAYERR                !" レイリー数の相対許容誤差
      REAL       RAYLL , RAYLU         !" レイリー数下限, 上限
*
      SAVE
*
      DATA       ITMAX         / 20 /
      DATA       RAYERR        / 1e-6 / 
      DATA       RAYLL, RAYLU  / 1e5, 1e6 /
*
      NAMELIST   /NMCRT/  ITMAX , RAYERR , RAYLU, RAYLL
*
      INTEGER    IFPAR, JFPAR
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
*
      IF ( OFIRST )THEN
         OFIRST = .FALSE.
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMCRT , END=1190 )
 1190    WRITE  ( JFPAR, NMCRT  )
      ENDIF
*
      IMAX = ITMAX
      ERR  = RAYERR
      RL   = RAYLL
      RU   = RAYLU
*
      RETURN
      END
**********************************************************************
* PACKAGE UBISEC  !" 分割法による方程式の解の推定
*
*"[HIS]  93/07/03(takepiro) Numerical Recipies Sec.9.3 より移植(ZBRENT) 
*
*"    3 点の関数値から方程式の解を推定する. 関数値を計算し本ルーチンを
*"    繰り返し適用することにより与えられた精度の解を決定する. 
*
**********************************************************************
      SUBROUTINE UBRENT              !" Brent 法による方程式の解の推定
     M              ( X1   , X2   , X   , 
     M                FX1  , FX2  , FX  , 
     O                OCNV ,
     I                TOL                  )
*
*" [MODIFY]
      REAL   X1           !" 解の範囲その 1
      REAL   X2           !" 解の範囲その 2
      REAL   X            !" 推定解
*
      REAL   FX1          !" X1 における関数の値
      REAL   FX2          !" X2 における関数の値
      REAL   FX           !" X  における関数の値
*
*" [OUTPUT]
      LOGICAL OCNV        !" 収束したか否か
*
*" [INPUT]
      REAL   TOL          !" 収束判定条件
*
*" [INTERNAL WORK]
      REAL   D, E, P, Q, R, S, XM, TOL1
      REAL   EPS          !" 浮動小数精度
      LOGICAL OFIRST
      DATA    OFIRST / .TRUE. /
*
      REAL    RR , RP
      REAL*8  DELTPQ
      REAL*8  DELTPD
      REAL    DELTPS
      DATA    DELTPQ / 1.0 D-28 /         !' ４倍精度の収束条件
      DATA    DELTPD / 1.0 D-14 /         !'   倍精度の収束条件
      DATA    DELTPS / 1.0 E-7  /         !'   単精度の収束条件
*
      SAVE
*
*"  < 0. 浮動小数点精度チェック >
*
      IF( OFIRST )THEN
         OFIRST = .FALSE.
         WRITE ( 6, * ) ' BRENT METHOD : 93/07/05'
*
         RR  = 1.
         RP = RR + DELTPQ
         IF ( RP .GT. RR ) THEN
            WRITE ( 6,* ) ' precision check .. ', RP
            WRITE ( 6,* ) ' REAL*16 is used in UBRENT'
            EPS = DELTPQ
         ELSE
            RR  = 1.
            RP = RR + DELTPD
            WRITE ( 6,* ) ' precision check .. ', RP
            IF ( RP .GT. RR ) THEN
               WRITE ( 6,* ) ' REAL*8 is used in UBRENT'
               EPS = DELTPD
            ELSE
               WRITE ( 6,* ) ' REAL*4 is used in UBRENT'
               EPS = DELTPS
            ENDIF
         ENDIF
*
      ENDIF
*
*"  < 1. 配置がえ >
*
      IF ( FX * FX2 .GT. 0. ) THEN
         X2  = X1
         FX2 = FX1
         D   = X - X1
         E   = D
      ENDIF
      IF(ABS(FX2).LT.ABS(FX)) THEN
         X1  = X
         X   = X2
         X2  = X1
         FX1 = FX
         FX  = FX2
         FX2 = FX1
      ENDIF
*
*"  < 2. 収束判定 >
*
      TOL1 = 2.  * EPS * ABS(X) + 0.5*TOL
      XM   = 0.5 * ( X2 - X )
*
      IF(  ABS(XM) .LE. TOL1 .OR. FX. EQ.0. )THEN
         OCNV = .TRUE.
         RETURN
      ENDIF
*
*"  < 3. 新しい解の推定 >
*
      IF(  ABS(E).GE.TOL1 .AND. ABS(FX1).GT.ABS(FX)  ) THEN
         S=  FX / FX1
         IF ( X1 .EQ. X2 ) THEN                  !" はさみうち法
            P = 2. * XM * S
            Q = 1.-S
         ELSE                                    !" 逆 2 次関数補間
            Q = FX1 / FX2
            R = FX  / FX2
            P = S *( 2.*XM*Q*(Q-R)-(X-X1)*(R-1.) )
            Q = (Q-1.)*(R-1.)*(S-1.)
         ENDIF
         IF(P.GT.0.) Q=-Q
         P=ABS(P)
         IF( 2.*P .LT. MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q)) ) THEN
            E = D                           !" 新しい推定値を受け入れる
            D = P / Q
         ELSE
            D = XM                          !" 2 分法による推定
            E = D
         ENDIF
      ELSE
         D=XM                               !" 2 分法による推定
         E=D
      ENDIF
*
      X1  = X                       !"  以前の推定値を区間の端点として保存
      FX1 = FX
*
      IF(  ABS(D) .GT. TOL1 ) THEN  !" 新推定値の計算
         X  = X + D
      ELSE
         X  = X + SIGN(TOL1,XM)
      ENDIF
*
      OCNV = .FALSE.
*
      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:5) .EQ. 'GATOR' )THEN
         ANORM = ABSMAX
     I         ( GATOR, IDIM, JDIM, KMAX, IMAX, JMAX )
      ELSE IF( HITEM(1:5) .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,JDIM
            DO 8100 I=1,IDIM
               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 = ABS(DATA( I,J,K ))
            ENDIF
 1100 CONTINUE
*
      RETURN
      END
