* 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] 





* PACKAGE  ZXT10L32 !" 格子点数，波数 X-Z
*
      INTEGER     IMAX           !" 経度格子数
      INTEGER     JMAX           !" 緯度格子数
      INTEGER     IDIM           !" 経度配列サイズ
      INTEGER     JDIM           !" 緯度配列サイズ
      INTEGER     MMAX           !" 最大東西波数
      INTEGER     LMAX           !" 最大南北波数
      INTEGER     NMAX           !" 最大全波数
      INTEGER     MINT           !" 東西波数間隔
      INTEGER     KMAX           !" 鉛直層数
      INTEGER     KDIM           !" 鉛直配列サイズ
      INTEGER     NMDIM          !" 水平波数成分数
      INTEGER     JMXHF          !" JMAX/2+1
      INTEGER     NWORK          !" ワーク領域大きさ
*
      PARAMETER ( IMAX  = 32,     JMAX =  1,
     &            IDIM  = IMAX+1, JDIM = JMAX,
     &            MMAX  = 10,     LMAX =  0,    NMAX = MMAX, MINT = 1,
     &            KMAX  = 32,     KDIM = KMAX,
     &            NMDIM = (MMAX/MINT+1)*(2*(NMAX+1)-MMAX)
     &                  - (NMAX-LMAX)/MINT*(NMAX-LMAX+1),
     &            JMXHF = JMAX/2+1,
     &            NWORK = IDIM*JDIM*KMAX*12                       )

* PACKAGE ZHDIM  !" 標準物理定数(common)  
*
*"  [HIS]  90/08/31(numaguti)
*
      ReaL*8       ER             !" 地球半径
      ReaL*8       GRAV                        !" 重力加速度
      ReaL*8       CP                          !" 大気定圧比熱
      ReaL*8       RAIR                        !" 大気気体定数
      ReaL*8       EL                          !" 水の凝結の潜熱
      ReaL*8       CPVAP                       !" 水蒸気定圧比熱
      ReaL*8       RVAP                        !" 水蒸気気体定数
      ReaL*8       DH2O                        !" 水の密度
      ReaL*8       EPSV                        !" 水と大気の分子量比ε
      ReaL*8       ES0                         !" ０℃飽和蒸気圧：Ｐａ
      ReaL*8       STB                         !" Stefan-Boltzman
      ReaL*8       FKARM                       !" Karman 定数
      ReaL*8       EPSVT                       !" 1/ε-1
*
      COMMON   /COMCON/
     &          ER    , GRAV  ,
     &          CP    , RAIR  ,
     &          EL    , CPVAP , RVAP  , DH2O  , EPSV  , ES0   ,
     &          STB   , FKARM , EPSVT

*
*   [OUTPUT] 
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) 
      ReaL*8       FLAPLA( NMDIM  )            
      ReaL*8       EDEL  ( NMDIM  )            
*"       : 球関数展開定数
*
      ReaL*8       CORIOL ( IDIM, JDIM )       
      ReaL*8       UVFACT ( IDIM, JDIM )       
*"       : 力学定数
*
      ReaL*8       DIFV  ( NMDIM  )            
      ReaL*8       DIFT  ( NMDIM  )            
*"       : 時間積分定数
*
      ReaL*8       ALAT  ( JDIM )              
      ReaL*8       DLAT  ( JDIM )              
*        : 座標格子定数
*
*   [INTERNAL WORK] 
      INTEGER    IDAYL                       
      INTEGER    IEFOLD                      
      INTEGER    IDORDR                      
*
      ReaL * 8   QSINLA( JDIM )              
      ReaL * 8   QGW   ( JDIM )              
      ReaL*8       PI                          
      ReaL*8       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 / DBLE( 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                              
     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] 




* PACKAGE  ZXT10L32 !" 格子点数，波数 X-Z
*
      INTEGER     IMAX           !" 経度格子数
      INTEGER     JMAX           !" 緯度格子数
      INTEGER     IDIM           !" 経度配列サイズ
      INTEGER     JDIM           !" 緯度配列サイズ
      INTEGER     MMAX           !" 最大東西波数
      INTEGER     LMAX           !" 最大南北波数
      INTEGER     NMAX           !" 最大全波数
      INTEGER     MINT           !" 東西波数間隔
      INTEGER     KMAX           !" 鉛直層数
      INTEGER     KDIM           !" 鉛直配列サイズ
      INTEGER     NMDIM          !" 水平波数成分数
      INTEGER     JMXHF          !" JMAX/2+1
      INTEGER     NWORK          !" ワーク領域大きさ
*
      PARAMETER ( IMAX  = 32,     JMAX =  1,
     &            IDIM  = IMAX+1, JDIM = JMAX,
     &            MMAX  = 10,     LMAX =  0,    NMAX = MMAX, MINT = 1,
     &            KMAX  = 32,     KDIM = KMAX,
     &            NMDIM = (MMAX/MINT+1)*(2*(NMAX+1)-MMAX)
     &                  - (NMAX-LMAX)/MINT*(NMAX-LMAX+1),
     &            JMXHF = JMAX/2+1,
     &            NWORK = IDIM*JDIM*KMAX*12                       )

*
*   [OUTPUT] 
      ReaL*8       FLAPLA( NMDIM  )            
      ReaL*8       EDEL  ( NMDIM  )            
*
*   [INPUT] 
      ReaL*8       ER                          
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) 
*
*   [INTERNAL WORK] 
      INTEGER    L, M, LEND
      ReaL*8       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] 




* PACKAGE  ZXT10L32 !" 格子点数，波数 X-Z
*
      INTEGER     IMAX           !" 経度格子数
      INTEGER     JMAX           !" 緯度格子数
      INTEGER     IDIM           !" 経度配列サイズ
      INTEGER     JDIM           !" 緯度配列サイズ
      INTEGER     MMAX           !" 最大東西波数
      INTEGER     LMAX           !" 最大南北波数
      INTEGER     NMAX           !" 最大全波数
      INTEGER     MINT           !" 東西波数間隔
      INTEGER     KMAX           !" 鉛直層数
      INTEGER     KDIM           !" 鉛直配列サイズ
      INTEGER     NMDIM          !" 水平波数成分数
      INTEGER     JMXHF          !" JMAX/2+1
      INTEGER     NWORK          !" ワーク領域大きさ
*
      PARAMETER ( IMAX  = 32,     JMAX =  1,
     &            IDIM  = IMAX+1, JDIM = JMAX,
     &            MMAX  = 10,     LMAX =  0,    NMAX = MMAX, MINT = 1,
     &            KMAX  = 32,     KDIM = KMAX,
     &            NMDIM = (MMAX/MINT+1)*(2*(NMAX+1)-MMAX)
     &                  - (NMAX-LMAX)/MINT*(NMAX-LMAX+1),
     &            JMXHF = JMAX/2+1,
     &            NWORK = IDIM*JDIM*KMAX*12                       )

*
*   [OUTPUT] 
      ReaL*8       DIFV  ( NMDIM  )            
      ReaL*8       DIFT  ( NMDIM  )            
*
*   [INPUT] 
      INTEGER    IEFOLD                      
      INTEGER    IDORDR                      
      INTEGER    NMO   ( 2, 0:MMAX, 0:LMAX ) 
*
*   [INTERNAL WORK] 
      ReaL*8       DIFMAX
      ReaL*8       DIFVAL
      INTEGER    L, M, NM, LEND
      ReaL*8       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. )/DBLE(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. /DBLE(NMAX) )**(IDORDR/2)
 1200    CONTINUE
      ENDIF
*"                                      
*
*"                 < 2. 拡散係数を決める >
*
      DIFVAL = 1. / DBLE(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
*
*"                                      
*
*"                 < 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
