* PACKAGE DINTG  
*
*"[履歴]  97/05/19 堀之内武 bug fix
***********************************************************************
      SUBROUTINE TINTGR         
     M         ( WDVOR , WDDIV , WDT   , WDPI  , WDQ   ,
     I           WTVOR , WTDIV , WTT   , WTPI  , WTQ   ,
     I           DMLU  , KPIV  , DELT  ,
     C           DIFV  , DIFT  , FLAPLA, WPHIS ,
     C           GG    , TMTX  , DSIG                   )
*
*"                 Arakawa & Suarez 鉛直差分 σ座標
*"                   semi implicit scheme
*"                   implicit horizontal diffusion
*




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

*
*   [MODIFY] 
      ReaL*8       WDVOR ( NMDIM , KMAX     )  
      ReaL*8       WDDIV ( NMDIM , KMAX     )  
      ReaL*8       WDT   ( NMDIM , KMAX     )  
      ReaL*8       WDPI  ( NMDIM            )  
      ReaL*8       WDQ   ( NMDIM , KMAX     )  
*
*   [INPUT] 
      ReaL*8       WTVOR ( NMDIM, KMAX      )  
      ReaL*8       WTDIV ( NMDIM, KMAX      )  
      ReaL*8       WTT   ( NMDIM, KMAX      )  
      ReaL*8       WTPI  ( NMDIM            )  
      ReaL*8       WTQ   ( NMDIM, KMAX      )  
*
      ReaL*8       DMLU  ( NMDIM, KMAX, KMAX ) 
      INTEGER    KPIV  ( NMDIM, KMAX       ) 
      ReaL*8       DELT                        
*
      ReaL*8       DIFV  ( NMDIM  )            
      ReaL*8       DIFT  ( NMDIM  )            
      ReaL*8       FLAPLA( NMDIM  )            
      ReaL*8       WPHIS ( NMDIM )             
*
      ReaL*8       GG    ( KMAX )              
      ReaL*8       TMTX  ( KMAX, KMAX )        
      ReaL*8       DSIG  ( KMAX )              
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           WDTX  , WDPIX ,
     &           WDPHIX, WBRDIV
      ReaL*8       WDTX     ( NMDIM, KMAX )    
      ReaL*8       WDPIX    ( NMDIM )          
      ReaL*8       WDPHIX   ( NMDIM, KMAX )    
      ReaL*8       WBRDIV   ( NMDIM, KMAX )    
*
      ReaL*8       DELT2
      INTEGER    NM, K, KK
*
      DELT2  = DELT*2.  
*
*"         < 1. 非重力波項の仮積分 >
*
      DO 1100 K = 1 , KMAX
         DO 1100 NM=1,NMDIM 
            WDTX  ( NM,K ) = WDT  ( NM,K )*( 1. - DELT*DIFT( NM ) )
     &                     + WTT  ( NM,K )* DELT
 1100 CONTINUE
*
      DO 1200 NM = 1 , NMDIM
            WDPIX ( NM )   = WDPI ( NM ) + WTPI  ( NM )* DELT
 1200 CONTINUE
*
*"         < 2. ジオポテンシャル >
*
      CALL WHYDRO
     O         ( WDPHIX,
     I           WDTX   )
*
*"         < 3. 発散方程式の右辺  >
*
      DO 3100 K = 1 , KMAX
*
         DO 3110 NM = 1 , NMDIM
*
            WBRDIV( NM,K )
     &       =  WDDIV( NM,K ) * ( 1. + DELT *DIFV( NM ) )
     &                        * ( 1. + DELT2*DIFT( NM ) )
     &           + DELT*(   WTDIV ( NM,K )
     &                    - FLAPLA( NM )
     &                      * (  WDPHIX( NM,K )
     &                           + ( 1. + DELT2*DIFT( NM ) )
     &                             * (  WPHIS ( NM )
     &                                + WDPIX ( NM )*GG( K ) )  )   )
*
 3110    CONTINUE
 3100 CONTINUE
*
*"         < 4. 時間平均のＤをLU行列で解く >
*
      CALL LUSOLV
     M         ( WBRDIV,
     I           DMLU  , KPIV  ,
     D           1     , NMDIM , KMAX   )
*
*"         < 5. 熱力, 地表気圧の時間変化項 >
*
      DO 5100 K = 1, KMAX
         DO 5100 KK = 1 , KMAX
            DO 5100 NM = 1, NMDIM
               WTT ( NM,K ) = WTT   ( NM,K )
     &                       - TMTX ( K,KK ) * WBRDIV( NM,KK )
 5100 CONTINUE
*
      DO 5200 KK = 1 , KMAX
         DO 5200 NM = 1, NMDIM
              WTPI ( NM   ) = WTPI  ( NM   )
     &                        - DSIG(   KK ) * WBRDIV( NM,KK )
 5200 CONTINUE
*
*"         < 7. 時間積分  拡散は前方差分 >
*
      DO 7100  K = 1 , KMAX
         DO 7100  NM = 1 , NMDIM
*
*"        ------------------  < 渦度 > ---------------------------
*
            WDVOR( NM,K ) = ( WDVOR( NM,K ) + WTVOR( NM,K )*DELT2 )
     &                     /( 1. + DELT2*DIFV( NM )               )
*
*"        ------------------  < 発散 >  -------------------------
*
            WDDIV ( NM,K ) =  WBRDIV( NM,K )*2. - WDDIV ( NM,K )
*
*"        ------------------  < 温度 >  -------------------------
*
            WDT  ( NM,K ) = ( WDT  ( NM,K ) + WTT  ( NM,K )*DELT2 )
     &                    / ( 1. + DELT2*DIFT( NM )               )
*
*"        ------------------  < 比湿 >  -------------------------
*
            WDQ  ( NM,K ) = ( WDQ  ( NM,K ) + WTQ  ( NM,K )*DELT2 )
     &                    / ( 1. + DELT2*DIFT( NM )               )
*
 7100 CONTINUE
*
      DO 7200 NM = 1 , NMDIM
*
*"        ------------------  < ln(地表気圧) >  --------------------
*
            WDPI ( NM   ) =   WDPI ( NM   ) + WTPI ( NM   )*DELT2
*
 7200 CONTINUE
*
      RETURN
      END
***********************************************************************
      SUBROUTINE MOVGD        
     O         ( GBU   , GBV   , GBT   , GBPS  , GBQ   ,
     O           GBVOR , GBDIV ,
     I           GAU   , GAV   , GAT   , GAPS  , GAQ   ,
     I           GAVOR , GADIV                          )
*





* 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       GBU   ( IDIM, JDIM, KMAX )  
      ReaL*8       GBV   ( IDIM, JDIM, KMAX )  
      ReaL*8       GBT   ( IDIM, JDIM, KMAX )  
      ReaL*8       GBPS  ( IDIM, JDIM       )  
      ReaL*8       GBQ   ( IDIM, JDIM, KMAX )  
      ReaL*8       GBVOR ( IDIM, JDIM, KMAX )  
      ReaL*8       GBDIV ( IDIM, JDIM, KMAX )  
*"       : 格子点データ(t-Δt) ->  格子点データ(t)
*
*   [INPUT] 
      ReaL*8       GAU   ( IDIM, JDIM, KMAX )  
      ReaL*8       GAV   ( IDIM, JDIM, KMAX )  
      ReaL*8       GAT   ( IDIM, JDIM, KMAX )  
      ReaL*8       GAPS  ( IDIM, JDIM       )  
      ReaL*8       GAQ   ( IDIM, JDIM, KMAX )  
      ReaL*8       GAVOR ( IDIM, JDIM, KMAX )  
      ReaL*8       GADIV ( IDIM, JDIM, KMAX )  
*"       : 格子点データ(t)
*
      CALL COPY ( GBU   , GAU   , IDIM*JDIM*KMAX )
      CALL COPY ( GBV   , GAV   , IDIM*JDIM*KMAX )
      CALL COPY ( GBT   , GAT   , IDIM*JDIM*KMAX )
      CALL COPY ( GBPS  , GAPS  , IDIM*JDIM      )
      CALL COPY ( GBQ   , GAQ   , IDIM*JDIM*KMAX )
      CALL COPY ( GBVOR , GAVOR , IDIM*JDIM*KMAX )
      CALL COPY ( GBDIV , GADIV , IDIM*JDIM*KMAX )
*
      RETURN
      END
