* PACKAGE BALOCLI  !" 初期値設定 （ 剛体回転傾圧 ）
*"                  各層で剛体回転する流れ ＋ 初期微小擾乱 を設定。
*
*"  [HIS]   93/09/18 （堀之内）
**********************************************************************
      SUBROUTINE INISET               !" 初期値設定
     O         ( GAU   , GAV   , GAT   , GAPS  , GAQ   ,
     I           DUM1  , DUM2  , ALAT  , DUM3  ,
     I           SIG   , DUM4  , SIGM  , DUM5  )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [OUTPUT] 
      REAL       GAU   ( IDIM, JDIM, KMAX ),
     &           GAV   ( IDIM, JDIM, KMAX ),
     &           GAT   ( IDIM, JDIM, KMAX ),
     &           GAPS  ( IDIM, JDIM       )
      REAL       GAQ   ( IDIM, JDIM, KMAX )
*
*   [INPUT]
      REAL       SIG   ( KMAX )           !" σレベル(整数)
      REAL       SIGM  (KMAX+1)           !" σレベル(半整数)
      REAL       ALAT  ( JDIM )           !" 緯度
      REAL       DUM1,DUM2,DUM3,DUM4,DUM5 !" ダミー

*   [INTERNAL PARM] 
      REAL       TB ( KMAX )              !" 基本温度(レベル平均)
      REAL       TB1,TB2                  !" 基本温度 at σ=1,σkmax
      DATA       TB1,TB2 / 280. , 220. /
      REAL       U0 ( KMAX )              !" 赤道での風速
      REAL       U01,U02                  !" 同上     at σ=1,σkmax
      DATA       U01,U02 / 0. , 50. /
      REAL       PS0                      !" 基本 Ps (hPa)
      DATA       PS0         / 1000. /
      REAL       SALPHA( KMAX )           !" 静水圧の式の係数 α
      REAL       SBETA ( KMAX )           !" 静水圧の式の係数 β
      REAL       OMEGA 
      SAVE
*
*   [INTERNAL WORK] 
      REAL       C ( JDIM )               !" COS( 緯度 )
      REAL       S ( JDIM )               !" SIN( 緯度 )
      REAL       TW( KMAX )               !" 温度の波成分の係数
*
      CALL PCONST        !" 標準物理定数(デフォルト) 
      OMEGA = 2 * 3.141593 / 86400.
*
      DO 10 J=1,JDIM
        C( J ) = COS( ALAT(J) )
        S( J ) = SIN( ALAT(J) )
   10 CONTINUE
*
      DO 20 K=1,KMAX
        U0( K ) = U01 + (U02-U01) * ( 1-SIG(K) ) / ( 1 - SIG(KMAX) )
        TB( K ) = TB1 + (TB2-TB1) * ( 1-SIG(K) ) / ( 1 - SIG(KMAX) )
   20 CONTINUE
*
      DO 100 K = 1,KMAX 
      DO 100 J = 1,JDIM
      DO 100 I = 1,IDIM
        GAU( I, J, K ) = U0( K ) * C( J ) 
        GAV( I, J, K ) = 0
  100 CONTINUE
*
      IF (  U01 .NE. 0 ) THEN
        WRITE(6,*) '(INISET) U01 MUST BE ZERO'
        STOP
      ENDIF
      DO 120 J = 1,JDIM
      DO 120 I = 1,IDIM
        GAPS( I, J ) = PS0
  120 CONTINUE
*
      CALL HYDROS         !" 静水圧の係数の設定(dynamics/dhydr.F)
     I         ( SIG   , SIGM  )
      CALL HYDROQ         !" 静水圧の係数の参照(dynamics/dhydr.F)
     O         ( SALPHA, SBETA )
*
      TW ( 1 ) = U0(1)/CP/SALPHA(1)
      DO 140 K = 2,KMAX
        TW ( K ) = ( U0(K)-U0(K-1) )/CP/SALPHA(K)
     &             - TW(K-1)*SBETA(K-1)/SALPHA(K)
  140 CONTINUE
*
      DO 160 K = 1,KMAX 
      DO 160 J = 1,JDIM
      DO 160 I = 1,IDIM
        GAT( I, J, K ) =  TB(K) + OMEGA*ER* ( 1./3 - S(J)**2 ) * TW(K)
  160 CONTINUE
*
      CALL JOURAN         !" 初期擾乱
     M       ( GAU, GAV )
*
      RETURN
      END
************************************************************************
*" 初期擾乱(Legendre各成分エネルギー一定, 位相ランダム／そろっている)
*
      SUBROUTINE JOURAN
     M       ( GAU, GAV )         !" 初期擾乱
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.F"                    !" 格子点数, 波数
#endif
*
*   [MODIFY] 
      REAL       GAU   ( IDIM, JDIM, KMAX ),
     &           GAV   ( IDIM, JDIM, KMAX )
*
*   [INTERNAL PARAM] 
      INTEGER    MMN,MMX, LMN,LMX
**      DATA       MMN,MMX, LMN,LMX/ 1,MMAX-4, 0,3 /
      REAL       UORDER   !" 各成分の SQRT(2*ENRGY) : 発散成分がないので
*"                        !  Potential Energyは０であることに注意
      DATA       UORDER/ 0.001 /
      SAVE
*
*   [INTERNAL WORK] 
      REAL       WDPSI( NMDIM , KMAX ) !" 流線関数 ψ/ER (m/s)
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
*
*   [CONSTS]
      MMN = 1
      MMX = MMAX-3
      LMN = 0
      LMX = LMAX 

      CALL LTINIT
      CALL LTGET
     I       ( 'NMO', 2*(MMAX+1)*(LMAX+1),
     O          NMO  )
*
      DO 100 K = 1, KMAX
      DO 100 M = MMN,MMX
      DO 100 L = LMN,LMX
        IF ( NMO(1,M,L).NE.0 .AND. NMO(2,M,L).NE.0 ) THEN
          AR = SIN( 1. * (L+1) * M ) !" 位相をバラバラにする。
          AI = COS( 1. * (L+1) * M ) !" 位相をバラバラにする。
**          AR = 1                     !" 位相をそろえる。
**          AI = 0                     !" 位相をそろえる。
          N  =  M + L
          C  =  1./SQRT( 1.* N*(N+1) )
          WDPSI( NMO(1,M,L) , K ) = AR * C * UORDER
          WDPSI( NMO(2,M,L) , K ) = AI * C * UORDER
        ENDIF
  100 CONTINUE
*
      CALL W2G
     O         ( GAU   ,
     I           WDPSI ,
     I           'YGRA', 'SUB ', KMAX )
*
      IF ( MMAX .GE. 1 ) THEN
         CALL W2G
     O         ( GAV   ,
     I           WDPSI ,
     I           'XGRA', 'ADD ', KMAX )
      ENDIF
*
      RETURN
      END
