* PACKAGE P2SFCM  !" 物理過程(2) 地表フラックス
*"                   Ishiwatari (97/10/31) 係数出力
***********************************************************************
      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 )   !" 温度  Ｔ
c      REAL       GDTM  ( IDIM*JDIM, KMAX )   !" 温度（半整数）
      REAL       GDTM  ( IDIM*JDIM, KMAX+1 )   !" 温度（半整数）
      REAL       GDTG  ( IDIM*JDIM       )   !" 地表面温度
      REAL       GDQ   ( IDIM*JDIM, KMAX )   !" 比湿  ｑ
      REAL       GDP   ( IDIM*JDIM, KMAX )   !" 気圧
c      REAL       GDPM  ( IDIM*JDIM, KMAX )   !" 気圧（半整数）
      REAL       GDPM  ( IDIM*JDIM, KMAX+1 )   !" 気圧（半整数）
      REAL       GDZ   ( IDIM*JDIM, KMAX )   !" 高度
*
      REAL       GRBET ( IDIM*JDIM )         !" 地表湿潤度β
      REAL       GRZ0  ( IDIM*JDIM )         !" 地表粗度長
      INTEGER    IDSFC ( IDIM*JDIM )         !" 地表状態
*
*   [INTERNAL PARM] 
      REAL       P00, THS
      REAL       URMIN
      REAL       USMINM, USMINH, USMINE
      REAL       USMAXM, USMAXH, USMAXE
      DATA       P00   / 100000. /
      DATA       THS   / 300.    /            !" 基本温位
      DATA       URMIN / 0.01    /            !" Ｒｉ数用風最小値
      DATA       USMINM/ 0.01    /            !" 運動量用風最小値
      DATA       USMINH/ 0.01    /            !" 熱用風最小値
      DATA       USMINE/ 0.01    /            !" 水蒸気用風最小値
      DATA       USMAXM/ 1000.   /            !" 運動量用風最大値
      DATA       USMAXH/ 1000.   /            !" 熱用風最大値
      DATA       USMAXE/ 1000.   /            !" 水蒸気用風最大値
*
      NAMELIST  /NMSFCM/ THS   , URMIN ,
     &                   USMINM, USMAXM,
     &                   USMINH, USMAXH,
     &                   USMINE, USMAXE
*
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           EXNER , VABS  , RIB   ,
     &           CMV   , CHV   , CEV   ,
     &           QSAT  , DQSAT
      REAL       EXNER ( IDIM*JDIM )         !" 大気Exner関数
      REAL       VABS  ( IDIM*JDIM )         !" 風速絶対値
      REAL       RIB   ( IDIM*JDIM )         !" バルクＲｉ数
      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 - P2.1a DATE=92/04/17'
         OFIRST = .FALSE.
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMSFCM, END=190 )
  190    WRITE  ( JFPAR, NMSFCM )
*
         CALL GZDBGT ( 'MISS', VMISS )
*
         CALL HISTRG
     I       ( 'CDM' , 'Bulk coef (mom)       ' ,'m^2/s    ', 'ASFC',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
         CALL HISTRG
     I       ( 'CDH' , 'Bulk coef (heat)      ' ,'m^2/s    ', 'ASFC',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
         CALL HISTRG
     I       ( 'CDE' , 'Bulk coef (q)         ' ,'m^2/s    ', 'ASFC',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
         CALL HISTRG
     I       ( 'CRM' , 'Bulk coef * rho( mom) ' ,'Kg m/s    ', 'ASFC',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
         CALL HISTRG
     I       ( 'CRH' , 'Bulk coef * rho (heat)' ,'Kg m/s    ', 'ASFC',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
         CALL HISTRG
     I       ( 'CRE' , 'Bulk coef * rho (q)   ' ,'Kg m/s    ', 'ASFC',
     I          VMISS  , VMISS , VMISS , VMISS , ISTYPL  ,
     I         'CON'   , '   ' , 0     , 0     , 'X'     ,'(1PE12.3)' )
*
      ENDIF
*
      AKAPPA = RAIR / CP
      DO 100 IJ = 1, IDIM*JDIM
         EXNER ( IJ ) = ( GDP ( IJ,1 )/GDPM( IJ,1 ) ) ** AKAPPA
  100 CONTINUE
*
*"         < 1. バルクＲｉ数 >
*
      DO 1100 IJ = 1, IDIM*JDIM
         VABS ( IJ ) = SQRT( GDU( IJ,1 )**2+GDV( IJ,1 )**2 )
         RIB  ( IJ ) = GRAV / THS
     &                * (   GDT ( IJ,1 )/EXNER ( IJ ) 
     &                    - GDTG( IJ   )              )
     &                / MAX( VABS( IJ ), URMIN )**2 * GDZ( IJ,1 )
 1100 CONTINUE
*
*"         < 3. バルク係数 >
*
      CALL BLKCOF
     O         ( CMV   , CHV   , CEV   ,
     I           RIB   , VABS  , GRZ0  , GRZ0  , GDZ   )
*
      CALL HISTIN ( CMV  , 'CDM'  )
      CALL HISTIN ( CHV  , 'CDH'  )
      CALL HISTIN ( CEV  , 'CDE'  )
*
*"         < 4. バルク輸送係数 >
*
      DO 4100 IJ = 1, IDIM*JDIM
         CMV  ( IJ ) = CMV ( IJ )
     &                 * GDPM( IJ,1 ) / ( RAIR * GDTM( IJ,1 ) )
     &                 * MIN( MAX ( VABS( IJ ), USMINM ), USMAXM )
         CHV  ( IJ ) = CHV ( IJ )
     &                 * GDPM( IJ,1 ) / ( RAIR * GDTM( IJ,1 ) )
     &                 * MIN( MAX ( VABS( IJ ), USMINH ), USMAXH )
         CEV  ( IJ ) = CEV ( IJ )
     &                 * GDPM( IJ,1 ) / ( RAIR * GDTM( IJ,1 ) )
     &                 * MIN( MAX ( VABS( IJ ), USMINE ), USMAXE )
 4100 CONTINUE
*
      CALL HISTIN ( CMV  , 'CRM'  )
      CALL HISTIN ( CHV  , 'CRH'  )
      CALL HISTIN ( CEV  , 'CRE'  )
*
*"         < 5. 飽和比湿  >
*
      DO 5100 IJ = 1, IDIM*JDIM
         QSAT ( IJ ) = FQSAT ( GDTG( IJ ), GDPM( IJ,1 ) )
         DQSAT( IJ ) = FDQSAT( GDTG( IJ ), QSAT( IJ )   )
 5100 CONTINUE
*
*"         < 6. フラックス >
*
      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 )
         TFLUX( IJ,1 ) = TFLUX( IJ,1 )
     &                   + CP *  CHV ( IJ )
     &                        *(   GDTG( IJ )
     &                           - GDT ( IJ,1 )/EXNER ( IJ )   )
         QFLUX( IJ,1 ) = QFLUX( IJ,1 )
     &                   + EL *  CEV ( IJ ) * GRBET( IJ )
     &                        *( QSAT( IJ ) - GDQ ( IJ,1 ) )
 6100 CONTINUE
*
*"       < 7. 陰解用行列 >
*
      CALL RESET ( TMTXS, IDIM*JDIM*6 )
      CALL RESET ( QMTXS, IDIM*JDIM*6 )
*
      DO 7100 IJ = 1, IDIM*JDIM
         UMTXS( IJ )      =         CMV( IJ )
*
         TMTXS( IJ,1,0  ) =    CP * CHV( IJ ) / EXNER( IJ )
         TMTXS( IJ,0,1  ) =  - CP * CHV( IJ ) / EXNER( IJ )
*
         QMTXS( IJ,1,0  ) =    CP * CEV( IJ ) * GRBET( IJ )
         QMTXS( IJ,0,1  ) =  - CP * CEV( IJ ) * GRBET( IJ )
*
         IF ( IDSFC( IJ ) .GE. 1 ) THEN
*
            TMTXS( IJ,1,-1 ) = - CP * CHV( IJ )
            TMTXS( IJ,0,0  ) =   CP * CHV( IJ )
*
            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
