* PACKAGE P2VDFM  !" 物理過程(2) 鉛直拡散
*
***********************************************************************
      SUBROUTINE VDFUVT        !" 鉛直拡散
     O         ( UFLUX , VFLUX , TFLUX , QFLUX ,
     M           UMTX  , TMTX  , QMTX  ,
     I           GDU   , GDV   , GDT   , GDTM  , GDQ   ,
     I           GDP   , GDPM  , GDZ   , GDZM           )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [OUTPUT] 
      REAL       UFLUX ( IDIM*JDIM, KMAX+1 ) !" Ｕのフラックス
      REAL       VFLUX ( IDIM*JDIM, KMAX+1 ) !" Ｖのフラックス
      REAL       TFLUX ( IDIM*JDIM, KMAX+1 ) !" Ｔのフラックス
      REAL       QFLUX ( IDIM*JDIM, KMAX+1 ) !" ｑのフラックス
*
*   [MODIFY] 
      REAL       UMTX  ( IDIM*JDIM ,   KMAX, -1:1 ) !" Ｕの陰解行列
      REAL       TMTX  ( IDIM*JDIM , 0:KMAX, -1:1 ) !" Ｔの陰解行列
      REAL       QMTX  ( IDIM*JDIM ,   KMAX, -1:1 ) !" ｑの陰解行列
*
*   [INPUT] 
      REAL       GDU   ( IDIM*JDIM, KMAX )   !" 西風  Ｕ
      REAL       GDV   ( IDIM*JDIM, KMAX )   !" 南風  Ｖ
      REAL       GDT   ( IDIM*JDIM, KMAX )   !" 温度  Ｔ
      REAL       GDTM  ( IDIM*JDIM, KMAX+1 ) !" 温度（半整数）
      REAL       GDQ   ( IDIM*JDIM, KMAX )   !" 比湿  ｑ
      REAL       GDP   ( IDIM*JDIM, KMAX )   !" 気圧  Ｐ
      REAL       GDPM  ( IDIM*JDIM, KMAX )   !" 気圧（半整数）
      REAL       GDZ   ( IDIM*JDIM, KMAX )   !" 高度
      REAL       GDZM  ( IDIM*JDIM, KMAX )   !" 高度（半整数）
*
*   [INTERNAL PARM] 
      REAL       P00, THS
      REAL       URMIN, RIBMIN
      DATA       P00   / 100000. /
      DATA       THS   / 300.    /           !" 基本温位
      DATA       URMIN / 0.1     /           !" 風二乗差最小値
      DATA       RIBMIN/ -100.   /           !" バルクＲｉ数最小値
*
      NAMELIST  /NMVDFM/ THS   , URMIN , RIBMIN
*
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE

*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           DVDZ  , RIB   ,
     &           DFH   , DFM   ,
     &           EXNER , EXNERM,
     &           WORK1 , WORK2 , WORK3 , WORK4
      REAL       DVDZ  ( IDIM*JDIM, KMAX )   !" d|v|/dz
      REAL       RIB   ( IDIM*JDIM, KMAX )   !" バルクＲｉ数
      REAL       DFH   ( IDIM*JDIM, KMAX+1 ) !" 輸送係数：Ｔ
      REAL       DFE   ( IDIM*JDIM, KMAX+1 ) !" 輸送係数：ｑ
      REAL       DFM   ( IDIM*JDIM, KMAX+1 ) !" 輸送係数：ｕ
      REAL       EXNER ( IDIM*JDIM, KMAX )   !" Exner関数：整数
      REAL       EXNERM( IDIM*JDIM, KMAX )   !" Exner関数：半整数
      REAL       WORK1 ( IDIM*JDIM, KMAX )   !" ワーク
      REAL       WORK2 ( IDIM*JDIM, KMAX )
      REAL       WORK3 ( IDIM*JDIM, KMAX )
      REAL       WORK4 ( IDIM*JDIM, KMAX )
*
      INTEGER    IJ, K
      INTEGER    IFPAR, JFPAR
      REAL       AKAPPA
      REAL       VMISS
*
*"         < 1. 定数 >
*
      IF ( OFIRST ) THEN
         WRITE ( 6,* ) ' VERTICAL DIFFUSION  P2.1 DATE=90/12/26'
         OFIRST = .FALSE.
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMVDFM, END=1190 )
 1190    WRITE  ( JFPAR, NMVDFM )
*
         CALL GZDBGT ( 'MISS', VMISS )
*
         CALL HISTRG
     I       ( 'KFM' , 'vert diff coef (mom) ' ,'m^2/s    ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
         CALL HISTRG
     I       ( 'KFH' , 'vert diff coef (heat) ' ,'m^2/s    ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
         CALL HISTRG
     I       ( 'KFE' , 'vert diff coef (q)    ' ,'m^2/s    ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
         CALL HISTRG
     I       ( 'DFM' , 'diff coef * rho (mom) ' ,'Kg m/s    ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
         CALL HISTRG
     I       ( 'DFH' , 'diff coef * rho (heat)' ,'Kg m/s    ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
         CALL HISTRG
     I       ( 'DFE' , 'diff coef * rho (q)   ' ,'Kg m/s    ', 'ALEV',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      ENDIF
*
      AKAPPA = RAIR / CP
      DO 1100 K = 1, KMAX
         DO 1110 IJ = 1, IDIM*JDIM
            EXNER ( IJ,K ) = ( GDP ( IJ,K )/P00 ) ** AKAPPA
            EXNERM( IJ,K ) = ( GDPM( IJ,K )/P00 ) ** AKAPPA
 1110    CONTINUE
 1100 CONTINUE
*
*"         < 2. バルクＲｉ数 >
*
      CALL RESET( DVDZ, IDIM*JDIM )
      CALL RESET( RIB , IDIM*JDIM )
*
      DO 2100 K = 2, KMAX
         DO 2110 IJ = 1, IDIM*JDIM
            DVDZ( IJ,K ) = SQRT( MAX( URMIN,
     &                               ( GDU( IJ,K )-GDU( IJ,K-1 ) )**2
     &                              +( GDV( IJ,K )-GDV( IJ,K-1 ) )**2
     &                         )    )
     &                       / ( GDZ( IJ,K ) - GDZ( IJ,K-1 )          )
*
            RIB ( IJ,K ) = GRAV / THS
     &                     * (   GDT( IJ,K   )/EXNER( IJ,K   )
     &                         - GDT( IJ,K-1 )/EXNER( IJ,K-1 )  )
     &                     / ( GDZ( IJ,K ) - GDZ( IJ,K-1 )      )
     &                     / DVDZ ( IJ,K )**2
*
            RIB ( IJ,K ) = MAX (  RIB( IJ,K ), RIBMIN  )
 2110    CONTINUE
 2100 CONTINUE
*
*"         < 3. 拡散係数 >
*
      CALL VDFCOF
     O         ( DFM   , DFH   , DFE   ,
     I           RIB   , DVDZ  , GDZM  ,
     W           WORK1 , WORK2 , WORK3 , WORK4 )
*
*"         < 3.1 浅い積雲対流 >
*
      CALL SHLCOF
     M         ( DFM   , DFH   , DFE   ,
     I           GDT   , GDQ   , GDP   , GDZ   ,
     W           WORK1 , WORK2 , WORK3 ,
     W           WORK4(1,1), WORK4(1,2) )
*
      CALL HISTIN ( DFM  , 'KFM'  )
      CALL HISTIN ( DFH  , 'KFH'  )
      CALL HISTIN ( DFE  , 'KFE'  )
*
*"         < 4. 輸送係数 >
*
      DO 4100 K = 2, KMAX
         DO 4110 IJ = 1, IDIM*JDIM
            DFM( IJ,K ) = DFM( IJ,K )
     &                    *  GDPM( IJ,K ) / RAIR / GDTM( IJ,K )
     &                    / ( GDZ( IJ,K ) - GDZ( IJ,K-1 ) )
*
            DFH( IJ,K ) = DFH( IJ,K )
     &                    *  GDPM( IJ,K ) / RAIR / GDTM( IJ,K )
     &                    / ( GDZ( IJ,K ) - GDZ( IJ,K-1 ) )
*
            DFE( IJ,K ) = DFE( IJ,K )
     &                    *  GDPM( IJ,K ) / RAIR / GDTM( IJ,K )
     &                    / ( GDZ( IJ,K ) - GDZ( IJ,K-1 ) )
*
 4110    CONTINUE
 4100 CONTINUE
*
      CALL HISTIN ( DFM  , 'DFM'  )
      CALL HISTIN ( DFH  , 'DFH'  )
      CALL HISTIN ( DFE  , 'DFE'  )
*
*"         < 5. フラックス >
*
      DO 5100 K = 2, KMAX
         DO 5110 IJ = 1, IDIM*JDIM
*
            UFLUX( IJ,K ) = UFLUX( IJ,K )
     &                     + DFM ( IJ,K )
     &                       * ( GDU ( IJ,K-1 ) - GDU ( IJ,K ) )
*
            VFLUX( IJ,K ) = VFLUX( IJ,K )
     &                      + DFM ( IJ,K )
     &                       * ( GDV ( IJ,K-1 ) - GDV ( IJ,K ) )
*
            TFLUX( IJ,K ) = TFLUX( IJ,K )
     &                      + CP * DFH ( IJ,K )  * EXNERM( IJ,K   )
     &                        * (  GDT( IJ,K-1 ) / EXNER ( IJ,K-1 )
     &                           - GDT( IJ,K   ) / EXNER ( IJ,K   ) )
*
            QFLUX( IJ,K ) = QFLUX( IJ,K )
     &                      + EL * DFE ( IJ,K )
     &                       * ( GDQ( IJ,K-1 ) - GDQ( IJ,K ) )
*
 5110    CONTINUE
 5100 CONTINUE
*
*"         < 6. 陰解用行列 >
*
      DO 6100 K = 2, KMAX
         DO 6110 IJ = 1, IDIM*JDIM
            UMTX( IJ,K,0  ) = UMTX( IJ,K, 0 ) + DFM  ( IJ,K )
            UMTX( IJ,K,-1 ) =                 - DFM  ( IJ,K )
*
            TMTX( IJ,K,0  ) = TMTX( IJ,K, 0 )
     &                          + CP * DFH ( IJ,K   )*EXNERM( IJ,K   )
     &                                               /EXNER ( IJ,K   )
            TMTX( IJ,K,-1 ) =   - CP * DFH ( IJ,K   )*EXNERM( IJ,K   )
     &                                               /EXNER ( IJ,K-1 )
*
            QMTX( IJ,K,0  ) = QMTX( IJ,K, 0 ) + CP * DFE  ( IJ,K )
            QMTX( IJ,K,-1 ) =                 - CP * DFE  ( IJ,K )
 6110    CONTINUE
 6100 CONTINUE
*
      DO 6200 K = 1, KMAX-1
         DO 6210 IJ = 1, IDIM*JDIM
            UMTX( IJ,K,0  ) = UMTX( IJ,K,0 ) + DFM ( IJ,K+1 )
            UMTX( IJ,K,1  ) =                - DFM ( IJ,K+1 )
*
            TMTX( IJ,K,0  ) = TMTX( IJ,K,0 )
     &                          + CP * DFH ( IJ,K+1 )*EXNERM( IJ,K+1 )
     &                                               /EXNER ( IJ,K   )
            TMTX( IJ,K,1  ) =   - CP * DFH ( IJ,K+1 )*EXNERM( IJ,K+1 )
     &                                               /EXNER ( IJ,K+1 )
*
            QMTX( IJ,K,0  ) = QMTX( IJ,K,0 ) + CP * DFE ( IJ,K+1 )
            QMTX( IJ,K,1  ) =                - CP * DFE ( IJ,K+1 )
 6210    CONTINUE
 6200 CONTINUE
*
      RETURN
      END
