* PACKAGE DSETC  !" 力学 定数設定(水平差分)
*
*"  [HIS] 90/08/31(numaguti)
*"        92/05/18(masaaki)     SETCOF: EDEL(LMAX=0)
*"        92/07/20(takepiro)    SETDIF              エラーチェック
*"        93/09/04(pochaka)     DYNPAR              設定順序変更
*" 
***********************************************************************
      SUBROUTINE DSETC     !"  力学パートの定数セット(水平差分)
     O         ( NMO   , FLAPLA, EDEL  ,
     O           CORIOL, UVFACT, DIFV  , DIFT  ,
     C           ALAT  , DLAT                   )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [OUTPUT] 
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       EDEL  ( NMDIM  )            !" ζ，Ｄ→Ｕ，Ｖ
*"       : 球関数展開定数
*
      REAL       CORIOL ( IDIM, JDIM )       !" コリオリ係数 ｆ
      REAL       UVFACT ( IDIM, JDIM )       !" u→U のファクター
*"       : 力学定数
*
      REAL       DIFV  ( NMDIM  )            !" 運動量水平拡散係数
      REAL       DIFT  ( NMDIM  )            !" 熱，水水平拡散係数
*"       : 時間積分定数
*
      REAL       ALAT  ( JDIM )              !" 緯度
      REAL       DLAT  ( JDIM )              !" 緯度荷重
*        : 座標格子定数
*
*   [INTERNAL WORK] 
      INTEGER    IDAYL                       !" 一日の長さ(秒)
      INTEGER    IEFOLD                      !" 最大波数減衰時間
      INTEGER    IDORDR                      !" 水平拡散のオーダー
*
      REAL * 8   QSINLA( JDIM )              !" sin（緯度）：倍精度
      REAL * 8   QGW   ( JDIM )              !" ガウス荷重 ：倍精度
      REAL       PI                          !" 円周率π
      REAL       OMEGA                       !" 自転角速度Ω
*
      INTEGER    I, J
*
*"         < 1. 力学パラメーター >
*
      CALL DYNPAR
     O         ( IDAYL ,
     O           IEFOLD, IDORDR )
*
*"         < 2. 球面調和関数展開パラメーター >
*
      IF ( JMAX .GE. 2 ) THEN
         CALL GAUSS
     O         ( QSINLA, QGW   ,
     I           JMAX            )
      ELSE
         QSINLA( 1 ) = SIN ( ALAT( 1 ) )
         QGW   ( 1 ) = DLAT( 1 )
      ENDIF
*
      CALL SPSTUP
     I         ( QSINLA, QGW   , ER    ,
     O           NMO                    )
*
      CALL SETCOF
     O         ( FLAPLA, EDEL  ,
     I           ER    , NMO     )
*
*"         < 3.  コリオリ因子, Uファクター )
*
      PI     = ATAN ( 1. )*4.
      OMEGA  = 2. *PI / REAL( IDAYL )
*
      DO 3100 I = 1, IDIM
         DO 3110 J = 1, JDIM
            CORIOL ( I,J ) = 2. *OMEGA * SIN( ALAT( J ) )
 3110    CONTINUE
 3100 CONTINUE
*
      DO 3200 I = 1, IDIM
         DO 3210 J = 1, JDIM
            UVFACT ( I,J ) = COS( ALAT( J ) )
 3210    CONTINUE
 3200 CONTINUE
*
*"         < 4. 水平拡散のファクター >
*
      CALL SETDIF
     O         ( DIFV  , DIFT  ,
     I           IEFOLD, IDORDR,
     I           NMO            )
*
      RETURN
      END
************************************************************************
      SUBROUTINE DYNPAR           !" 力学パラメーター
     O         ( IDAYL ,
     O           IEFOLD, ORDERY  )
*
*   [OUTPUT] 
      INTEGER    IDAYL                       !" 一日の長さ(秒)
      INTEGER    IEFOLD                      !" 最大波数減衰時間
      INTEGER    ORDERY                      !" 水平拡散のオーダー
*
*   [INTERNAL WORK] 
      INTEGER    DAYLEN                      !" 一日の長さ
      INTEGER    TEFOLD                      !" 最大波数減衰時間
      INTEGER    ORDER                       !" 水平拡散のオーダー
      CHARACTER  TUNIT*4                     !" 時間の単位
*
      INTEGER    IFPAR, JFPAR
*
      NAMELIST    /NMDAYL/ TUNIT, DAYLEN
      NAMELIST    /NMHDIF/ TUNIT, TEFOLD, ORDER
*
*"  <  一日の長さの設定. >
*
      TUNIT  = 'HOUR'
      DAYLEN = 24
*
      CALL   REWNML ( IFPAR , JFPAR )
      READ   ( IFPAR, NMDAYL, END=1190 )
 1190 WRITE  ( JFPAR, NMDAYL )
*
      CALL GUCT2I                              !" 保坂 93/09/04 移動
     I         ( DAYLEN, TUNIT  , 'SEC' , 1     ,
     O           IDAYL                                  )
*
*"  <  最大波数減衰時間の設定 >
*
      TUNIT  = 'HOUR'
      TEFOLD = 4
      ORDER  = 16
*
      CALL   REWNML ( IFPAR , JFPAR )
      READ   ( IFPAR, NMHDIF, END=1290 )
 1290 WRITE  ( JFPAR, NMHDIF )
*
      CALL GUCT2I
     I         ( TEFOLD, TUNIT  , 'SEC' , 1     ,
     O           IEFOLD                                 )
*
*"  <  水平拡散の次数の設定 > 
*
      ORDERY = ORDER
*
      RETURN
      END
***********************************************************************
      SUBROUTINE SETCOF              !" スペクトル計算の係数
     O        ( FLAPLA, EDEL  ,
     I          ER    , NMO   )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT] 
      REAL       FLAPLA( NMDIM  )            !" ラプラシアンの係数
      REAL       EDEL  ( NMDIM  )            !" ζ,Ｄ→Ｕ,Ｖ
*
*   [INPUT] 
      REAL       ER                          !" 地球半径
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
*
*   [INTERNAL WORK] 
      INTEGER    L, M, LEND
      REAL       AN
*
*"         < 1. ラプラシアンとその変形 >
*
      DO 2110 M = 0 , MMAX, MINT
         LEND = MIN( LMAX, NMAX-M )
         DO 2100 L = 0 , LEND
            AN = L + M
            IF ( AN .GT. 0.1 ) THEN
               FLAPLA( NMO(1,M,L) )= - AN*( AN+1. ) / ER**2
               FLAPLA( NMO(2,M,L) )= - AN*( AN+1. ) / ER**2
               EDEL  ( NMO(1,M,L) )= - ER / ( AN *( AN+1. ) )
               EDEL  ( NMO(2,M,L) )= - ER / ( AN *( AN+1. ) )
            ENDIF
 2100    CONTINUE
 2110 CONTINUE
*
      FLAPLA( NMO(1,0,0) ) =  0.  
      FLAPLA( NMO(2,0,0) ) =  0.  
      IF ( LMAX .NE. 0 ) THEN
         EDEL  ( NMO(1,0,0) ) =  0.  
         EDEL  ( NMO(2,0,0) ) =  0.  
      ELSE
         EDEL  ( NMO(1,0,0) ) =  ER
         EDEL  ( NMO(2,0,0) ) =  0.  
      ENDIF
*
      RETURN
      END
***********************************************************************
      SUBROUTINE SETDIF               !" 水平拡散のファクター
     O         ( DIFV  , DIFT  ,
     I           IEFOLD, IDORDR,
     I           NMO             )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#endif
*
*   [OUTPUT] 
      REAL       DIFV  ( NMDIM  )            !" 運動量水平拡散係数
      REAL       DIFT  ( NMDIM  )            !" 熱，水水平拡散係数
*
*   [INPUT] 
      INTEGER    IEFOLD                      !" 最大波数減衰時間
      INTEGER    IDORDR                      !" 水平拡散のオーダー
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) !" スペクトルの添字順番
*
*   [INTERNAL WORK] 
      REAL       DIFMAX
      REAL       DIFVAL
      INTEGER    L, M, NM, LEND
      REAL       AN, DIFTN
*
*"                 < 1. 拡散係数をかけない部分 >
*
      DIFMAX = 0.  
*
      DO 1100 M = 0 , MMAX, MINT
         LEND = MIN( LMAX, NMAX-M )
         DO 1110 L = 0 , LEND
            AN = L + M
            IF ( AN .GT. 0.1 ) THEN
               DIFTN  = ( AN*( AN+1. )/REAL(NMAX) )**(IDORDR/2)
               DIFMAX = MAX( DIFMAX, DIFTN )
               DIFT( NMO(1,M,L) )  = DIFTN
               DIFT( NMO(2,M,L) )  = DIFTN
            ENDIF
 1110    CONTINUE
 1100 CONTINUE
*
      IF( NMAX .GT. 0 )THEN
         DO 1200 NM = 1 , NMDIM
            DIFV( NM ) = DIFT( NM ) - ( 2. /REAL(NMAX) )**(IDORDR/2)
 1200    CONTINUE
      ENDIF
*"                                      !" 竹広 92/07/20 上記 IF 節追加
*
*"                 < 2. 拡散係数を決める >
*
      DIFVAL = 1. / REAL(IEFOLD)
*
*"                 < 3. 拡散のファクター >
*
      IF( DIFMAX .GT. 0.0 )THEN           
         DO 3100 NM = 1 , NMDIM
            DIFV ( NM ) = DIFVAL * DIFV( NM ) / DIFMAX
            DIFT ( NM ) = DIFVAL * DIFT( NM ) / DIFMAX
 3100    CONTINUE
      ELSE
         CALL MSGDMP('W','SETDIF',' HORIZONTAL DIFF. COEF. IS ZERO. ' )
      ENDIF
*
*"                                      !" 竹広 92/07/20 上記 IF 節追加
*
*"                 < 4. 拡散しない成分 >
*
      DIFT ( NMO( 1,0,0 ) ) = 0.  
      DIFT ( NMO( 2,0,0 ) ) = 0.  
      DIFV ( NMO( 1,0,0 ) ) = 0.  
      DIFV ( NMO( 2,0,0 ) ) = 0.  
      IF ( LMAX. GE. 1 ) THEN
         DIFV ( NMO( 1,0,1 ) ) = 0.  
         DIFV ( NMO( 2,0,1 ) ) = 0.  
      ENDIF
      IF ( ( MMAX. GE. 1 ).AND.( MINT .EQ. 1 ) ) THEN
         DIFV ( NMO( 1,1,0 ) ) = 0.  
         DIFV ( NMO( 2,1,0 ) ) = 0.  
      ENDIF
*
      RETURN
      END
