* PACKAGE REST2   !" 初期値設定 （ 静止、温度場は米標準大気1976 
*"                                 + 極微小温度擾乱）
*"                       ： 単精度のとき擾乱は消えてしまうことに注意
*
*"  [HIS]   93/11/05 （堀之内）
**********************************************************************
      SUBROUTINE INISET               !" 初期値設定
     O         ( GAU   , GAV   , GAT   , GAPS  , GAQ   ,
     I           DUM1  , DUM2  , DUM3  , DUM4  ,
     I           SIG   , DUM6  , DUM7  , DUM8  )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数, 波数
#else
#include        "zcdim.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 )              !" σレベル(整数)
*
*   [INTERNAL PARAM] 
      SAVE
      REAL       TPRTRB                      !" 温度擾乱振幅
      DATA       TPRTRB / 1.D-10 /
      INTEGER    NNODE
      PARAMETER (NNODE = 18 )
      REAL       PNODE (NNODE)               !" 気圧 : 0 - 120 km
      DATA       PNODE/1013.   , 540.5   , 265.0   , 121.1   , 55.29 , 
     &                 25.49   , 11.97   , 5.746   , 2.871   , 1.491 ,
     &                 0.7908  , 0.2196  , 0.05221 , 0.01052 , 1.836E-3,
     &                 3.201E-4, 7.104E-5, 2.538E-5 /
      REAL       TNODE (NNODE)               !" 気温 : 0 - 120 km
      DATA       TNODE/288.15  , 255.68  , 223.25  , 216.65  , 216.65  ,
     &                 221.55  , 226.51  , 236.51  , 250.35  , 264.16  ,
     &                 270.65  , 247.02  , 219.59  , 198.64  , 186.97  ,
     &                 195.08  , 240.00  , 360.00   /
*
*   [INTERNAL WORK] 
      REAL       RANDOM( IDIM, JDIM, KMAX )
      INTEGER    I, J, K, KN
      REAL       GTB (KMAX)
      REAL       Z, ZN1, ZN2
*
*
      CALL RESET ( GAU   , IDIM*JDIM*KMAX )
      CALL RESET ( GAV   , IDIM*JDIM*KMAX )
      CALL RESET ( GAT   , IDIM*JDIM*KMAX )
      CALL RESET ( GAQ   , IDIM*JDIM*KMAX )
      CALL RESET ( GAPS  , IDIM*JDIM      )
*
*" << INIT : これが初期値だ >> 
*
      DO 1000 K = 1, KMAX
      DO 1000 J = 1, JDIM
      DO 1000 I = 1, IDIM
         GAU  ( I,J,K ) = 0.
         GAV  ( I,J,K ) = 0.
         GAQ  ( I,J,K ) = 0.
 1000 CONTINUE
*
      DO 1500 K = 1, KMAX
        Z = -LOG( SIG(K) )
        DO 1510 KN = 1,NNODE
          ZN1 = -LOG( PNODE(KN)  /PNODE(1) )
          ZN2 = -LOG( PNODE(KN+1)  /PNODE(1) )
          IF(Z.GE.ZN1 .AND. Z.LT.ZN2) THEN
            GTB(K) = RINTRP(Z, ZN1,ZN2, TNODE(KN), TNODE(KN+1) )
            GOTO 1520
          ENDIF
 1510   CONTINUE
        IF(Z.GE.ZN2) GTB(K)=RINTRP(Z,ZN2,ZN2,TNODE(NNODE),TNODE(NNODE))
        IF(Z.LT.0.) CALL MSGDMP('E','INISET','soitsu ha hen dayo')
 1520   CONTINUE
        DO 1600 J = 1, JDIM
        DO 1600 I = 1, IDIM
          GAT  ( I,J,K ) = GTB(K)
 1600   CONTINUE
 1500 CONTINUE
*
      DO 2000 J = 1, JDIM
      DO 2000 I = 1, IDIM
         GAPS ( I,J   ) = 1000.
 2000 CONTINUE
*
*" << RANDOM : 乱数の付加 >> 
*
      DO 3000 K = 1, KMAX
      DO 3000 J = 1, JDIM
      DO 3000 I = 1, IDIM
         RANDOM( I,J,K )
     &        =  SIN(  REAL(I*J*(KMAX-K))
     &               / REAL(IMAX*JMAX*KMAX)*10. )
 3000 CONTINUE
*
      DO 3100 K = 1, KMAX
      DO 3100 J = 1, JDIM
      DO 3100 I = 1, IDIM
         GAT ( I,J,K ) = GAT ( I,J,K )
     &                   + RANDOM( I,J,K ) * TPRTRB
 3100 CONTINUE
*
      RETURN
      END
**********************************************************************
      REAL FUNCTION RINTRP(Z, Z1,Z2, R1,R2 )
      REAL      Z,Z1,Z2,R1,R2, T
*
      T = (Z-Z1)/(Z2-Z1)
      RINTRP = (1-T)*R1 + T*R2
*
      RETURN
      END
