* PACKAGE P2SFCM  !" 物理過程(2) 地表フラックス(拡散係数一定)
*
*  [HIS] 92/07/20(takepiro)
*
***********************************************************************
      SUBROUTINE SFCUVT           !"  地表フラックス
     M         ( UFLUX , VFLUX , TFLUX , QFLUX ,
     M           UMTXS , TMTXS , QMTXS ,
     I           GDU   , GDV   , GDT   , GDTM  , GDTG  , GDQ   ,
     I           GDP   , GDPM  , GDZ   ,
     I           GRBET , GRZ0  , IDSFC                          )
*
*   [PARAM] 
#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 ) !" ｑのフラックス
*
      REAL       UMTXS ( IDIM*JDIM )             !" ｕ陰解行列：地表
      REAL       TMTXS ( IDIM*JDIM , 0:1, -1:1 ) !" Ｔ陰解行列：地表
      REAL       QMTXS ( IDIM*JDIM , 0:1, -1:1 ) !" ｑ陰解行列：地表
*
*   [INPUT] 
      REAL       GDU   ( IDIM*JDIM, KMAX )   !" 西風  ｕ
      REAL       GDV   ( IDIM*JDIM, KMAX )   !" 南風  ｖ
      REAL       GDT   ( IDIM*JDIM, KMAX )   !" 温度  Ｔ
      REAL       GDTM  ( IDIM*JDIM, KMAX )   !" 温度（半整数）
      REAL       GDTG  ( IDIM*JDIM       )   !" 地表面温度
      REAL       GDQ   ( IDIM*JDIM, KMAX )   !" 比湿  ｑ
      REAL       GDP   ( IDIM*JDIM, KMAX )   !" 気圧
      REAL       GDPM  ( IDIM*JDIM, KMAX )   !" 気圧（半整数）
      REAL       GDZ   ( IDIM*JDIM, KMAX )   !" 高度
*
      REAL       GRBET ( IDIM*JDIM )         !" 地表湿潤度β
      REAL       GRZ0  ( IDIM*JDIM )         !" 地表粗度長
      INTEGER    IDSFC ( IDIM*JDIM )         !" 地表状態
*
*   [INTERNAL PARM] 
      REAL       P00
      DATA       P00   / 100000. /
*
      REAL       DFMCON, DFHCON, DFECON
      DATA       DFMCON /    0.0   /         !" ｕ拡散係数デフォルト
      DATA       DFHCON /    0.0   /         !" Ｔ拡散係数デフォルト
      DATA       DFECON /    0.0   /         !" ｑ拡散係数デフォルト
*
      NAMELIST  /NMVDF0/ DFMCON, DFHCON, DFECON
*
      LOGICAL    OTFFIX                      !" 温度フラックス固定 SW
      DATA       OTFFIX / .FALSE. /
      REAL       TGFLUX                      !" 地表面熱フラックス
      DATA       TGFLUX /    0.0   /         !" TGFLUX デフォルト
*
      NAMELIST   /NMTFLX/ OTFFIX, TGFLUX
*
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           EXNER , 
     &           CMV   , CHV   , CEV   ,
     &           QSAT  , DQSAT
      REAL       EXNER ( IDIM*JDIM )         !" 大気Exner関数
      REAL       CMV   ( IDIM*JDIM )         !" バルク輸送係数：ｕ
      REAL       CHV   ( IDIM*JDIM )         !" バルク輸送係数：Ｔ
      REAL       CEV   ( IDIM*JDIM )         !" バルク輸送係数：ｑ
      REAL       QSAT  ( IDIM*JDIM )         !" 地表飽和比湿
      REAL       DQSAT ( IDIM*JDIM )         !" 地表飽和比湿変化
*
      INTEGER    IJ
      INTEGER    IFPAR, JFPAR
      REAL       AKAPPA
*
#ifdef SYS_IBMS
      INCLUDE   (ZQSAT)                      !" 飽和比湿
#else
#include        "zqsat.F"                    !" 飽和比湿
#endif
*
      IF ( OFIRST ) THEN
         WRITE ( 6,* ) ' SURFACE FLUX - CONST.DIFF. DATE=92/07/14'
         OFIRST = .FALSE.
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMVDF0, END=190 )
  190    WRITE  ( JFPAR, NMVDF0 )
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMTFLX, END=290 )
  290    WRITE  ( JFPAR, NMTFLX )
*
      ENDIF
*
      AKAPPA = RAIR / CP
      DO 100 IJ = 1, IDIM*JDIM
         EXNER ( IJ ) = ( GDP ( IJ,1 )/GDPM( IJ,1 ) ) ** AKAPPA
  100 CONTINUE
*
*"         < 1. 拡散係数 >
*
      CALL SETV( CMV, DFMCON, IDIM*JDIM )
      CALL SETV( CHV, DFHCON, IDIM*JDIM )
      CALL SETV( CEV, DFECON, IDIM*JDIM )
*
*"         < 2. 輸送係数 >
*
      DO 4100 IJ = 1, IDIM*JDIM
         CMV  ( IJ ) = CMV ( IJ )
     &                 * GDPM( IJ,1 ) / ( RAIR * GDTM( IJ,1 ) )
     &                    / ( GDZ( IJ,1 ) - 0.0 )
         CHV  ( IJ ) = CHV ( IJ )
     &                 * GDPM( IJ,1 ) / ( RAIR * GDTM( IJ,1 ) )
     &                    / ( GDZ( IJ,1 ) - 0.0 )
         CEV  ( IJ ) = CEV ( IJ )
     &                 * GDPM( IJ,1 ) / ( RAIR * GDTM( IJ,1 ) )
     &                    / ( GDZ( IJ,1 ) - 0.0 )
 4100 CONTINUE
*
*"         < 3. 飽和比湿  >
*
      DO 5100 IJ = 1, IDIM*JDIM
         QSAT ( IJ ) = FQSAT ( GDTG( IJ ), GDPM( IJ,1 ) )
         DQSAT( IJ ) = FDQSAT( GDTG( IJ ), QSAT( IJ )   )
 5100 CONTINUE
*
*"         < 4. フラックス >
*
      DO 6100 IJ = 1, IDIM*JDIM
         UFLUX( IJ,1 ) = UFLUX( IJ,1 )
     &                         - CMV ( IJ ) * GDU( IJ,1 )
         VFLUX( IJ,1 ) = VFLUX( IJ,1 )
     &                         - CMV ( IJ ) * GDV( IJ,1 )
         IF ( OTFFIX )THEN
            TFLUX( IJ,1 ) = TFLUX( IJ,1 ) + TGFLUX
         ELSE
            TFLUX( IJ,1 ) = TFLUX( IJ,1 )
     &                        + CP *  CHV ( IJ )
     &                        *(   GDTG( IJ )
     &                           - GDT ( IJ,1 )/EXNER ( IJ )   )
         ENDIF
         QFLUX( IJ,1 ) = QFLUX( IJ,1 )
     &                   + EL *  CEV ( IJ ) * GRBET( IJ )
     &                        *( QSAT( IJ ) - GDQ ( IJ,1 ) )
 6100 CONTINUE
*
*"       < 5. 陰解用行列 >
*
      CALL RESET ( TMTXS, IDIM*JDIM*6 )
      CALL RESET ( QMTXS, IDIM*JDIM*6 )
*
      DO 7100 IJ = 1, IDIM*JDIM
         UMTXS( IJ )      =         CMV( IJ )
*
         IF( OTFFIX )THEN
            TMTXS( IJ,1,0  ) = 0.
            TMTXS( IJ,0,1  ) = 0.
         ELSE
            TMTXS( IJ,1,0  ) =    CP * CHV( IJ ) / EXNER( IJ )
            TMTXS( IJ,0,1  ) =  - CP * CHV( IJ ) / EXNER( IJ )
         ENDIF
*
         QMTXS( IJ,1,0  ) =    CP * CEV( IJ ) * GRBET( IJ )
         QMTXS( IJ,0,1  ) =  - CP * CEV( IJ ) * GRBET( IJ )
*
         IF ( IDSFC( IJ ) .GE. 1 ) THEN
*
            IF( OTFFIX )THEN
               TMTXS( IJ,1,-1 ) = 0.
               TMTXS( IJ,0,0  ) = 0.
            ELSE
               TMTXS( IJ,1,-1 ) = - CP * CHV( IJ )
               TMTXS( IJ,0,0  ) =   CP * CHV( IJ )
            ENDIF
*
            QMTXS( IJ,1,-1 ) = - EL * CEV( IJ ) * GRBET( IJ )
     &                                          * DQSAT( IJ )
            QMTXS( IJ,0,0  ) =   EL * CEV( IJ ) * GRBET( IJ )
     &                                          * DQSAT( IJ )
         ENDIF
*
 7100 CONTINUE
*
      RETURN
      END
