* PACKAGE P2IMTX  !" 物理過程(2) 陰解行列処理
*
***********************************************************************
      SUBROUTINE PHMTX0             !" 陰解用の配列初期化
     O         ( UMTX  , TMTX  , QMTX  ,
     O           UFLUX , VFLUX , TFLUX , QFLUX ,
     I           GDPM  , DELT  , GRHC  , IDSFC  )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [OUTPUT] 
      REAL       UMTX  ( IDIM*JDIM,   KMAX, -1:1 ) !" Ｕの陰解行列
      REAL       TMTX  ( IDIM*JDIM, 0:KMAX, -1:1 ) !" Ｔの陰解行列
      REAL       QMTX  ( IDIM*JDIM,   KMAX, -1:1 ) !" ｑの陰解行列
      REAL       UFLUX ( IDIM*JDIM, KMAX+1 )     !" Ｕのフラックス
      REAL       VFLUX ( IDIM*JDIM, KMAX+1 )     !" Ｖのフラックス
      REAL       TFLUX ( IDIM*JDIM, KMAX+1 )     !" Ｔのフラックス
      REAL       QFLUX ( IDIM*JDIM, KMAX+1 )     !" ｑのフラックス
*
*   [INPUT] 
      REAL       GDPM  ( IDIM*JDIM, KMAX+1 ) !" 気圧（半整数）
      REAL       DELT                        !" 時間刻みΔt
      REAL       GRHC  ( IDIM*JDIM )         !" 地表熱容量
      INTEGER    IDSFC ( IDIM*JDIM )         !" 地表状態
*
*   [INTERNAL WORK] 
      INTEGER    IJ, K
*
*
*"         < 1. 質量, 熱容量の項 >
*
      CALL RESET ( UMTX, IDIM*JDIM*KMAX    *3  )
      CALL RESET ( TMTX, IDIM*JDIM*(KMAX+1)*3  )
      CALL RESET ( QMTX, IDIM*JDIM*KMAX    *3  )
*
*
      DO 1100 K = 1, KMAX
         DO 1100 IJ = 1, IDIM*JDIM
            UMTX ( IJ,K,0 ) =  ( GDPM ( IJ,K ) - GDPM( IJ,K+1 ) )
     &                          / GRAV / DELT
            TMTX ( IJ,K,0 ) =  UMTX ( IJ,K,0 ) * CP
            QMTX ( IJ,K,0 ) =  UMTX ( IJ,K,0 ) * CP
 1100 CONTINUE
*
      DO 1200 IJ = 1, IDIM*JDIM
         IF ( IDSFC( IJ ) .GE. 1 ) THEN
            TMTX ( IJ,0,0 ) = GRHC( IJ ) / DELT
         ELSE
            TMTX ( IJ,0,0 ) = 1.   
         ENDIF
 1200 CONTINUE
*
*"         < 2. フラックスをリセット >
*
      CALL RESET( UFLUX  , IDIM*JDIM*(KMAX+1) )
      CALL RESET( VFLUX  , IDIM*JDIM*(KMAX+1) )
      CALL RESET( TFLUX  , IDIM*JDIM*(KMAX+1) )
      CALL RESET( QFLUX  , IDIM*JDIM*(KMAX+1) )
*
      RETURN
      END
************************************************************************
      SUBROUTINE PHYTND   !" 時間変化率の計算 (implicit)
     O         ( GTU   , GTV   , GTT   , GTTG  , GTQ   ,
     I           UFLUX , VFLUX ,
     I           TFLUX , SFLUX , RFLUX , GFLUX , QFLUX ,
     I           UMTX  , TMTX  , QMTX  ,
     I           UMTXS , TMTXS , QMTXS ,
     I           RMTXS , DELT  , IDSFC                  )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [OUTPUT] 
      REAL       GTU   ( IDIM*JDIM, KMAX )   !" 東西運動量変化項ＵＡ
      REAL       GTV   ( IDIM*JDIM, KMAX )   !" 南北運動量変化項ＶＡ
      REAL       GTT   ( IDIM*JDIM, KMAX )   !" 温度時間変化項  Ｈ
      REAL       GTTG  ( IDIM*JDIM )         !" 地表温度変化率
      REAL       GTQ   ( IDIM*JDIM, KMAX )   !" 比湿時間変化項  Ｒ
*
*   [INPUT] 
      REAL       UFLUX ( IDIM*JDIM, KMAX+1 ) !" Ｕのフラックス
      REAL       VFLUX ( IDIM*JDIM, KMAX+1 ) !" Ｖのフラックス
      REAL       TFLUX ( IDIM*JDIM, KMAX+1 ) !" Ｔのフラックス
      REAL       SFLUX ( IDIM*JDIM  )        !" 地表短波フラックス
      REAL       RFLUX ( IDIM*JDIM  )        !" 地表長波フラックス
      REAL       GFLUX ( IDIM*JDIM  )        !" 地中フラックス
      REAL       QFLUX ( IDIM*JDIM, KMAX+1 ) !" ｑのフラックス
*
      REAL       UMTX  ( IDIM*JDIM,   KMAX, -1:1 ) !" ｕ陰解行列
      REAL       TMTX  ( IDIM*JDIM, 0:KMAX, -1:1 ) !" Ｔ陰解行列
      REAL       QMTX  ( IDIM*JDIM,   KMAX, -1:1 ) !" ｑ陰解行列
      REAL       UMTXS ( IDIM*JDIM            ) !" ｕ陰解行列：地表
      REAL       TMTXS ( IDIM*JDIM, 0:1, -1:1 ) !" Ｔ陰解行列：地表
      REAL       QMTXS ( IDIM*JDIM, 0:1, -1:1 ) !" ｑ陰解行列：地表
      REAL       RMTXS ( IDIM*JDIM,      -1:1 ) !" Ｔ陰解行列：放射
*
      REAL       DELT                        !" 時間刻みΔt
*
      INTEGER    IDSFC ( IDIM*JDIM )         !" 地表状態
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           GTTQX , ATQMTX
      REAL       GTTQX ( IDIM*JDIM, -KMAX:KMAX )      !" Ｔｑ時間変化
      REAL       ATQMTX( IDIM*JDIM, -KMAX:KMAX, -1:1 )!" ＬＵ行列
      REAL       AUVMTX( IDIM*JDIM,       KMAX, -1:1 )!" ＬＵ行列
      EQUIVALENCE  ( ATQMTX, AUVMTX )
*
      INTEGER    IJ, K, L
*
*"         < 1. U, V の解 >
*
      CALL COPY ( AUVMTX, UMTX, IDIM*JDIM*KMAX*3 )
*
      DO 1100 IJ = 1, IDIM*JDIM
         AUVMTX( IJ,1,0 ) = AUVMTX( IJ,1,0 ) + UMTXS( IJ )
 1100 CONTINUE
*
      CALL LUMAK3
     O         ( AUVMTX    ,
     I           IDIM*JDIM , KMAX  )
*
      DO 1200 K = 1, KMAX
         DO 1200 IJ = 1, IDIM*JDIM
            GTU   ( IJ, K ) = UFLUX ( IJ,K ) - UFLUX ( IJ,K+1 )
            GTV   ( IJ, K ) = VFLUX ( IJ,K ) - VFLUX ( IJ,K+1 )
 1200 CONTINUE
*
      CALL LUSOL3
     M         ( GTU   ,
     I           AUVMTX,
     D           1     , IDIM*JDIM, KMAX )
*
      CALL LUSOL3
     M         ( GTV   ,
     I           AUVMTX,
     D           1     , IDIM*JDIM, KMAX )
*
*"         < 2. T,q の解 >
*
      DO 2100 L = -1, 1
         DO 2110 K = 1, KMAX
            DO 2110 IJ = 1, IDIM*JDIM
               ATQMTX( IJ, K, L ) = TMTX ( IJ,K,L )
               ATQMTX( IJ,-K,-L ) = QMTX ( IJ,K,L )
 2110    CONTINUE
         DO 2120 IJ = 1, IDIM*JDIM
               ATQMTX( IJ, 1, L ) = TMTX ( IJ,1,L ) + TMTXS( IJ,1,L )
               ATQMTX( IJ,-1,-L ) = QMTX ( IJ,1,L ) + QMTXS( IJ,1,L )
 2120    CONTINUE
 2100 CONTINUE
*
      DO 2200 IJ = 1, IDIM*JDIM
               ATQMTX( IJ, 0, 0 ) = TMTX ( IJ,0,0 )
     &                            + TMTXS( IJ,0,0 ) + QMTXS( IJ,0,0 )
     &                            + RMTXS( IJ,  0 )
               ATQMTX( IJ, 0, 1 ) = TMTXS( IJ,0,1 )
     &                            + RMTXS( IJ,  1 )
               ATQMTX( IJ, 0,-1 ) = QMTXS( IJ,0,1 )
 2200 CONTINUE
*
      CALL LUMAK3
     O         ( ATQMTX    ,
     I           IDIM*JDIM , 2*KMAX+1  )
*
      DO 2300 K = 1, KMAX
         DO 2300 IJ = 1, IDIM*JDIM
            GTTQX ( IJ, K ) = TFLUX ( IJ,K ) - TFLUX ( IJ,K+1 )
            GTTQX ( IJ,-K ) = QFLUX ( IJ,K ) - QFLUX ( IJ,K+1 )
 2300 CONTINUE
*
      DO 2400 IJ = 1, IDIM*JDIM
         GTTQX ( IJ,0 ) = - SFLUX ( IJ   ) - RFLUX ( IJ   )
     &                    - TFLUX ( IJ,1 ) - QFLUX ( IJ,1 )
     &                    + GFLUX ( IJ   )
 2400 CONTINUE
*
      CALL LUSOL3
     M         ( GTTQX ,
     I           ATQMTX,
     D           1     , IDIM*JDIM, 2*KMAX+1 )
*
*"         < 3. 時間変化率 >
*
      DO 3100 K = 1, KMAX
         DO 3100 IJ = 1, IDIM*JDIM
            GTU   ( IJ,K ) = GTU   ( IJ, K ) / DELT
            GTV   ( IJ,K ) = GTV   ( IJ, K ) / DELT
            GTT   ( IJ,K ) = GTTQX ( IJ, K ) / DELT
            GTQ   ( IJ,K ) = GTTQX ( IJ,-K ) / DELT / EL * CP
 3100 CONTINUE
*
      DO 3200 IJ = 1, IDIM*JDIM
         IF ( IDSFC( IJ ) .GE. 1 ) THEN
            GTTG  ( IJ   ) = GTTQX ( IJ, 0 ) / DELT
         ELSE
            GTTG  ( IJ   ) = 0.  
         ENDIF
 3200 CONTINUE
*
      RETURN
      END
************************************************************************
      SUBROUTINE FLXCOR        !" フラックスの補正
     M         ( UFLUX , VFLUX , TFLUX , QFLUX ,
     I           GTU   , GTV   , GTT   , GTTG  , GTQ   ,
     I           UMTX  , TMTX  , QMTX  ,
     I           UMTXS , TMTXS , QMTXS ,
     I           DELT                                     )
*
*   [INPUT] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [MODIFY] 
      REAL       UFLUX ( IDIM*JDIM, KMAX+1 ) !" Ｕのフラックス
      REAL       VFLUX ( IDIM*JDIM, KMAX+1 ) !" Ｖのフラックス
      REAL       TFLUX ( IDIM*JDIM, KMAX+1 ) !" Ｔのフラックス
      REAL       QFLUX ( IDIM*JDIM, KMAX+1 ) !" ｑのフラックス
*
*   [INPUT] 
      REAL       GTU   ( IDIM*JDIM, KMAX )   !" 東西運動量変化項ＵＡ
      REAL       GTV   ( IDIM*JDIM, KMAX )   !" 南北運動量変化項ＶＡ
      REAL       GTT   ( IDIM*JDIM, KMAX )   !" 温度時間変化項  Ｈ
      REAL       GTTG  ( IDIM*JDIM )         !" 地表温度変化率
      REAL       GTQ   ( IDIM*JDIM, KMAX )   !" 比湿時間変化項  Ｒ
*
      REAL       UMTX  ( IDIM*JDIM,   KMAX, -1:1 ) !" ｕ陰解行列
      REAL       TMTX  ( IDIM*JDIM, 0:KMAX, -1:1 ) !" Ｔ陰解行列
      REAL       QMTX  ( IDIM*JDIM,   KMAX, -1:1 ) !" ｑ陰解行列
      REAL       UMTXS ( IDIM*JDIM            ) !" ｕ陰解行列：地表
      REAL       TMTXS ( IDIM*JDIM, 0:1, -1:1 ) !" Ｔ陰解行列：地表
      REAL       QMTXS ( IDIM*JDIM, 0:1, -1:1 ) !" ｑ陰解行列：地表
      REAL       DELT                           !" 時間刻みΔt
*
*   [INTERNAL WORK] 
      INTEGER    IJ, K
      REAL       ELF
*
      ELF = EL/CP
*
      DO 1100 K = 2, KMAX
         DO 1100 IJ = 1, IDIM*JDIM
            UFLUX ( IJ,K ) = UFLUX  ( IJ,K )
     &                     - (  UMTX( IJ,K  ,-1)* GTU( IJ,K-1 )
     &                        - UMTX( IJ,K-1, 1)* GTU( IJ,K   ) ) * DELT
*
            VFLUX ( IJ,K ) = VFLUX  ( IJ,K )
     &                     - (  UMTX( IJ,K  ,-1)* GTV( IJ,K-1 )
     &                        - UMTX( IJ,K-1, 1)* GTV( IJ,K   ) ) * DELT
*
            TFLUX ( IJ,K ) = TFLUX  ( IJ,K )
     &                     - (  TMTX( IJ,K  ,-1)* GTT( IJ,K-1 )
     &                        - TMTX( IJ,K-1, 1)* GTT( IJ,K   ) ) * DELT
*
            QFLUX ( IJ,K ) = QFLUX  ( IJ,K )
     &                     - (  QMTX( IJ,K  ,-1)* GTQ( IJ,K-1 )
     &                        - QMTX( IJ,K-1, 1)* GTQ( IJ,K   ) ) * DELT
     &                                                            * ELF
 1100 CONTINUE
*
      DO 2200 IJ = 1, IDIM*JDIM
            UFLUX ( IJ,1 ) = UFLUX   ( IJ,1 )
     &                        - UMTXS( IJ )* GTU( IJ,1 ) * DELT
*
            VFLUX ( IJ,1 ) = VFLUX   ( IJ,1 )
     &                       -  UMTXS( IJ )* GTV( IJ,1 ) * DELT
*
            TFLUX ( IJ,1 ) = TFLUX   ( IJ,1 )
     &                     - (  TMTXS( IJ,1,-1 ) * GTTG( IJ )
     &                        + TMTXS( IJ,1,0  ) * GTT ( IJ,1 ) ) * DELT
*
            QFLUX ( IJ,1 ) = QFLUX   ( IJ,1 )
     &                     - (  QMTXS( IJ,1,-1 ) * GTTG( IJ )
     &                        + QMTXS( IJ,1,0  ) * GTQ ( IJ,1 )
     &                                           * ELF          ) * DELT
 2200 CONTINUE
*
      RETURN
      END
