* PACKAGE ASETS  !" 共通 定数(ｒ座標値)
*
*"  [HIS] 92/09/17(takepiro)
*"        93/04/30(takepiro) 軸ファイル名称変更
*"        93/05/27(takepiro) 重み計算
*"        93/07/27(takepiro) 軸ファイル名称変更
*
***********************************************************************
      SUBROUTINE SETRAD      !"  ｒレベル 
     O         ( ARAD  , DRAD  , HARAD ,
     D           KDIM                    )
*
*"               (* チェビシェフポイント ｒ座標 *)
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCCOM)                      !" 標準物理定数
      INCLUDE   (ZHDIM)                      !" 文字列文字数
#else
#include        "zccom.F"                    !" 標準物理定数
#include        "zhdim.F"                    !" 文字列文字数
#endif
      INTEGER    KDIM
*
*   [OUTPUT] 
      REAL       ARAD  ( 0:KDIM )              !" ｒレベル
      REAL       DRAD  ( 0:KDIM )              !" Δｒ
      CHARACTER  HARAD  *(*)                 !" 鉛直軸名称
*
*    [INTERNAL WORK]    
      INTEGER    KDMAX
      PARAMETER  ( KDMAX = 100 )
      REAL       XSI  ( 0:KDMAX )              !" ξレベル
      REAL       DXSI ( 0:KDMAX )              !" Δξ
      INTEGER    K
*
      WRITE ( 6,* ) ' CHEBYCHEV RADIAL LEVEL : DATE=93/07/29'
*
*"         < 1. 名称セット >
*
      CALL STNRAD
     O         ( HARAD ,
     I           KDIM   )
*
*"         < 2. ｒ, Δｒの計算 >
*
      CALL SETXSI
     O         ( XSI , DXSI , 
     I           KDIM         )
*
      DO 3100 K = 0 , KDIM
        ARAD  ( K ) = ( ( ERI + ERO ) + XSI(K) * ( ERO - ERI ) )*0.5
 3100 CONTINUE
*
      DO 3200 K = 0, KDIM
        DRAD  ( K ) = DXSI( K ) * ( ERO - ERI ) * 0.5 * ARAD(K)**2
 3200 CONTINUE
*
      RETURN
      END
***********************************************************************
      SUBROUTINE SETXSI     !" チェビシェフポイントの計算
     O              (  XSI , DXSI , 
     I                 KDIM         )
*
*   [PARAM] 
      INTEGER    KDIM
*
*   [OUTPUT] 
      REAL       XSI  ( 0:KDIM )               !" ξレベル
      REAL       DXSI ( 0:KDIM )               !" ξレベル
*
*    [INTERNAL WORK]    
      REAL       PI                          !" 円周率
      INTEGER    J, K
*
*"  < 1. チェビシェフメッシュポイントの計算 > 
*
      PI = ATAN( 1.0D0 ) * 4.0D0
*
      DO 1000 K = 0, KDIM
         XSI ( K ) = COS ( PI * K / KDIM )
 1000 CONTINUE 
*
*"  < 2. 重みの計算 > 
*
      DO 2100 K = 0, KDIM
         DXSI( K ) = 0. 
         DO 2200 J = 0, KDIM, 2
            DXSI( K ) = DXSI( K ) 
     &                  + 2. / (1-J**2) * COS ( J*K*PI / KDIM )
 2200    CONTINUE 
         DXSI( K ) = DXSI( K ) 
     &               - 1 
     &              -  0.5 * ( (-1)**(KDIM-1) -1 ) / ( 1 - KDIM**2 )
         DXSI( K ) = 2. / KDIM * DXSI( K )
 2100 CONTINUE 
*
      DXSI( 0 )    =  DXSI( 0 ) / 2 
      DXSI( KDIM ) =  DXSI( KDIM ) / 2 
*
      RETURN
      END
***********************************************************************
      SUBROUTINE STNRAD
     O               ( HARAD ,
     I                 KDIM   )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZHDIM)                      !" 文字列文字数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zhdim.F"                    !" 文字列文字数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
      INTEGER    KDIM
*
*   [OUTPUT] 
      CHARACTER  HARAD  *(*)             !" 鉛直軸名称：整数
*
*    [INTERNAL WORK]    
      INTEGER    INUM, IAZ
      CHARACTER  HNUM   *(NCC)
*
      CHARACTER  HARAD1 *(NCC)
      DATA       HARAD1 / 'CRAD' /       !" Chebychev Radial coordinate
*
      LOGICAL    LEPSL
*
*    [EXTERNAL FUNC]
      INTEGER    LENC
      LOGICAL    LREQ
*
*"         < 1. 名称セット >
*
      CALL GULCHR
     I     ( '(I4)', KDIM  ,
     O       HNUM  , INUM   )
*
      IAZ    = LENC( HARAD1 )
      HARAD  = HARAD1(1:IAZ)//HNUM(1:INUM)
*
      CALL GLPGET( 'LEPSL', LEPSL )
      IF ( .NOT.LEPSL )THEN
         CALL GLPSET( 'LEPSL', .TRUE. )
      ENDIF
*
      IF( LREQ( ERO-ERI,1.0 ) )THEN
         WRITE ( HNUM, '(F6.2)' ) ETA
         CALL CLADJ ( HNUM )
         INUM = LENC( HNUM )
         IAZ    = LENC( HARAD )
         HARAD  = HARAD(1:IAZ)//'-'//HNUM(1:INUM)
      ELSE IF ( LREQ( ERO,1.0 ) )THEN
         WRITE ( HNUM, '(F6.2)' ) ETA
         CALL CLADJ ( HNUM )
         INUM = LENC( HNUM )
         IAZ    = LENC( HARAD )
         HARAD  = HARAD(1:IAZ)//'N-'//HNUM(1:INUM)
      ELSE
         WRITE ( HNUM, '(F6.2)' ) ERI
         CALL CLADJ ( HNUM )
         INUM = LENC( HNUM )
         IAZ    = LENC( HARAD )
         HARAD  = HARAD(1:IAZ)//'R'//HNUM(1:INUM)
         WRITE ( HNUM, '(F6.2)' ) ERO
         CALL CLADJ ( HNUM )
         INUM = LENC( HNUM )
         IAZ    = LENC( HARAD )
         HARAD  = HARAD(1:IAZ)//'-'//HNUM(1:INUM)
      ENDIF
*
      IF ( .NOT.LEPSL )THEN
         CALL GLPSET( 'LEPSL', .FALSE. )
      ENDIF
*
      RETURN
      END