!= Module HeatFlux
!
! Authors::   ODAKA Masatsugu, TAKAHASHI Yoshiyuki
! Version::   $Id: surfaceflux_bulk.f90,v 1.1 2012-10-01 06:25:03 odakker Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006-2012. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!

module Surfaceflux_bulk
  !
  != 下部境界でのフラックスの計算モジュール
  !
  != Surface flux
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! Louis et al. (1982) の方法に基づいて地表面フラックスを計算. 
  !
  ! Surface fluxes are calculated by using the scheme by Louis et al. (1982). 
  !
  !== References
  !
  ! Louis, J-F., M. Tiedtke, and J-F. Geleyn, 
  ! A short history of the PBL parameterization at ECMWF, 
  ! Workshop on Planetary Boundary Layer Parameterization, 59-80, ECMWF, Reading, U.K., 
  ! 1982.
  !== Variable List
  !
  !== Procedures List
  !

  ! モジュール引用; USE statement
  !  

  ! GTOOL 変数と手続き
  ! GTOOL variables and procedures
  !
  use dc_types, only: DP, STRING
  use dc_iounit, only: FileOpen
  use dc_message, only: MessageNotify
  use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut

  ! 並列処理用変数
  ! Parallel processing variable 
  !       
  Use mpi_wrapper,only: myrank

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imin, & ! x 方向配列下限 
    &                      & ! Upper limit of array in x
    &                imax, & ! x 方向配列上限 
    &                      & ! Lower limit of array in x
    &                jmin, & ! y 方向配列下限 
    &                      & ! Upper limit of array in y
    &                jmax, & ! y 方向配列上限 
    &                      & ! Lower limit of array in y
    &                kmin, & ! z 方向配列下限 
    &                      & ! Upper limit of array in z
    &                kmax, & ! z 方向配列上限 
    &                      & ! Lower limit of array in z
    &                nx,   & ! x 方向格子点数 
    &                      & ! Number of grid point in x
    &                ny,   & ! y 方向格子点数 
    &                      & ! Number of grid point in y
    &                nz,   & ! z 方向格子点数 
    &                      & ! Number of grid point in z
    &                ncmax   ! 組成の数       
                             ! Number of spices

  use gridset_surfaceflux, only: GridsetSurfacefluxInit

  ! 座標軸と演算子
  ! Axes and operator settings
  !
  use axesset, only: z_dz,        & ! z 方向格子点間隔 
    &                             & ! Grid size in z
    &                xyz_Z,       & ! z 座標
    &                             & ! z coordinate
    &                xyz_avr_pyz, & ! 平均操作
    &                             & ! Average operator
    &                xyz_avr_xqz, & ! 平均操作
    &                             & ! Average operator
    &                xyr_avr_xyz, & ! 平均操作
    &                             & ! Average operator
    &                pyz_avr_xyz, & ! 平均操作
    &                             & ! Average operator
    &                xqz_avr_xyz, & ! 平均操作
    &                             & ! Average operator
    &                pyr_avr_xyr, & ! 平均操作
    &                             & ! Average operator
    &                xqr_avr_xyr    ! 平均操作
                                    ! Average operator

  ! 基本場変数
  ! Basic state variables
  !
  use basicset, only: xyz_ExnerBZ,  & ! 圧力関数
    &                               & ! Exner function
    &                 xyz_PressBZ,  & ! 圧力
    &                               & ! Pressure
    &                 xyz_PTempBZ,  & ! 温位
    &                               & ! Potential temperature
    &                 xyz_TempBZ,   & ! 温度
    &                               & ! Temperature
    &                 xyzf_QMixBZ,  & ! 混合比
    &                               & ! Mixing ration
    &                 xyz_DensBZ      ! 密度    
                                      ! Density

  ! 定数
  ! Constatns
  !
  use constants, only: Grav,       & ! 重力加速度
    &                              & ! Gravity
    &                  MolWtDry,   & ! 乾燥空気分子量 
    &                              & ! Molecular weight of dry air
    &                  PressBasis, & ! 温位の基準圧力
    &                              & ! Reference pressure
    &                  TempSfc,    & ! 地表面温度
    &                              & ! Surface temperature
    &                  PressSfc,   & ! 地表面気圧
    &                              & ! Surface pressure
    &                  CpDry,      & ! 乾燥空気定圧比熱 
    &                              & ! Specific heat of dry air
    &                  GasRDry       ! 乾燥空気気体定数
                                     ! Gas constant of dry air

  use constants0, only: FKarm        ! カルマン定数
                                     ! Karmann constant

  ! 組成に関する変数
  ! Setting for atmospheric composition
  ! 
  use composition, only: IdxCG,       & ! 凝結過程(気体)の配列添字 
    &                                 & ! Index of vapor
    &                    IdxCC,       & ! 凝結過程(雲)の配列添字 
    &                                 & ! Index of cloud
    &                    IdxG,        & ! 凝結過程(蒸気)の配列添字 
    &                                 & ! Index of vapor
    &                    SpcWetID,    & !凝結成分の化学種ID
    &                                 & ! ID number of moist condensation spices
    &                    CondNum,     & ! 雲の数
    &                                 & ! Number of cloud component
    &                    MolWtWet,    & ! 凝結成分分子量
    &                                 & ! Molecular weight of moist air
    &                    SpcWetSymbol   ! 凝結成分の化学種名
                                        ! Chemical formula of condensation spices

  ! 化学量計算に関する設定
  ! Settings for chemical calculation
  !
  use chemcalc, only: SvapPress ! 飽和蒸気圧
                                ! Saturation vapor pressure

  ! NAMELIST に関する設定
  ! Settings for NAMELIST
  !
  use namelist_util, only: namelist_filename ! NAMELIST ファイル名
                                             ! NAMELIST file name
  ! 時刻に関する設定
  ! Setting for time
  !
  use timeset, only:  TimeN ! 時刻 t 
                            ! Time "t"

  ! 時間変化項
  ! Tendency of variable
  use DExnerDt, only: xy_DExnerDt_xy_xyf ! D$\pi$/Dt

  ! dcpam 関連モジュール
  use surfaceflux_bulk_core, only: SurfaceFluxInit_Core => SurfaceFluxInit, &
                                    SurfaceFlux_Core => SurfaceFlux

  ! 暗黙の型宣言禁止
  ! Implicit none
  !
  implicit none

  ! 属性の指定
  !
  private

  ! 公開手続き
  ! Public procedure
  !
  public surfaceflux_bulk_init
  public surfaceflux_bulk_forcing

  ! 非公開変数
  ! Privete variables
  ! 
  logical, save:: FlagConstBulkCoef
                            ! Flag for using constant bulk coefficient
  logical, save:: FlagUseOfBulkCoefInNeutralCond
                            ! Flag for using bulk coefficient in neutral condition
  real(DP), save:: ConstBulkCoef
                            ! バルク係数一定値. 
                            ! Steady value of bulk coefficient
  real(DP), save :: VelMinForRi = 1.0d-8 
                            ! リチャード数計算用速度下限値
                            ! Lower limit of velocity for Ri
  real(DP), save :: SfcRoughLength = 1.0d-2
                            ! 祖度長さ
                            ! Roughness length
  real(DP), save :: VelBulkCoefMin = 0.0d0
                            ! $ u $ バルク係数最小値. 
                            ! Minimum value of $ u $ bulk coefficient
  real(DP), save :: TempBulkCoefMin = 0.0d0
                            ! $ T $ バルク係数最小値. 
                            ! Minimum value of $ T $ bulk coefficient
  real(DP), save :: QmixBulkCoefMin = 0.0d0
                            ! $ q $ バルク係数最小値. 
                            ! Minimum value of $ q $ bulk coefficient
  real(DP), save :: VelBulkCoefMax = 1.0d2
                            ! $ u $ バルク係数最大値. 
                            ! Maximum value of $ u $ bulk coefficient
  real(DP), save :: TempBulkCoefMax = 1.0d2
                            ! $ T $ バルク係数最大値. 
                            ! Maximum value of $ T $ bulk coefficient
  real(DP), save :: QmixBulkCoefMax = 1.0d2
                            ! $ q $ バルク係数最大値. 
                            ! Maximum value of $ q $ bulk coefficient

  real(DP), save :: PTempFluxMin = 0.0d0
                            ! 温位フラックス最小値. 
                            ! Minimum value of potential temp. flux
  real(DP), save :: ExnerFluxMin = 0.0d0
                            ! 圧力関数フラックス最小値. 
                            ! Minimum value of exner function flux
  real(DP), save :: QmixFluxMin = 0.0d0
                            ! 混合比フラックス最小値. 
                            ! Minimum value of mixing ratio flux
  real(DP), save  :: Vel0 = 0.0d0  ! 下層での水平速度嵩上げ値
                                   ! 

  character(*), parameter:: module_name = 'surfaceflux_bulk'
                                   ! モジュールの名称.
                                   ! Module name

contains
!!!------------------------------------------------------------------------!!!
  subroutine Surfaceflux_Bulk_init
    !
    ! NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. 
    !

    ! 暗黙の型宣言禁止
    ! Implicit none
    !
    implicit none

    ! 作業変数
    ! Work variables
    !
    integer    :: l,   &! 組成方向に回る DO ループ用作業変数
      &                & ! Work variables for DO loop in dimension of constituents
      &           unit   ! 出力装置番号
                         ! Device number 

    !---------------------------------------------------------------
    ! 格子点の初期化
    !
    call GridsetSurfacefluxInit

    !---------------------------------------------------------------
    ! NAMELIST から情報を取得
    !
!    NAMELIST /surface_flux_bulk_nml/ &
!      &  SfcRoughLength, Vel0                 

!    call FileOpen(unit, file=namelist_filename, mode='r')
!    read(unit, NML=surface_flux_bulk_nml)
!    close(unit)  

    ! dcpam ソースの初期化プログラムを call
    call SurfaceFluxInit_Core

!    if (myrank == 0) then 
!      call MessageNotify( "M", module_name, "SfcRoughLength = %f",  &
!        &  d=(/SfcRoughLength/))
!      call MessageNotify( "M", module_name, "Vel0 = %f", d=(/Vel0/))
!    end if

    call HistoryAutoAddVariable(      &
      & varname='PTempSfcFlux',       &
      & dims=(/'x','y','t'/),         &
      & longname='surface potential temperature flux (heat flux divided by density and specific heat)', &
      & units='K.m.s-1',             &
      & xtype='float')

    call HistoryAutoAddVariable(      &
      & varname='ExnerSfcFlux',       &
      & dims=(/'x','y','t'/),         &
      & longname='surface exner function flux (heat flux divided by density and specific heat)', &
      & units='s-1',             &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelXSfcFlux',    &
      & dims=(/'x','y','t'/),     &
      & longname='surface flux of x-component of velocity (momentum flux divided by density)', &
      & units='m2.s-2',           &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelYSfcFlux',    &
      & dims=(/'x','y','t'/),     &
      & longname='surface flux of y-component of velocity (momentum flux divided by density)', &
      & units='m2.s-2',           &
      & xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_SfcFlux', & 
        & dims=(/'x','y','t'/),     &
        & longname='surface flux of '          &
        &           //trim(SpcWetSymbol(l))//' mixing ratio (mass flux divided by density)',  &
        & units='m.s-1',    &
        & xtype='float')
    end do


    call HistoryAutoAddVariable(  &
      & varname='PTempSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='potential temperature tendency by surface flux', &
      & units='K.s-1',            &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='ExnerSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='exner function tendency by surface flux', &
      & units='K.s-1',            &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelXSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='x-component velocity tendency by surface flux', &
      & units='m.s-2',            &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='VelYSfc',         &
      & dims=(/'x','y','z','t'/), &
      & longname='y-component velocity tendency by surface flux', &
      & units='m.s-2',            &
      & xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_Sfc', & 
        & dims=(/'x','y','z','t'/),     &
        & longname=trim(SpcWetSymbol(l))//' mixing ratio tendency by surface flux',  &
        & units='s-1',    &
        & xtype='float')
    end do


    call HistoryAutoAddVariable(       &
      & varname='SfcHeatFlux',         &
      & dims=(/'x','y','t'/),          &
      & longname='surface heat flux',  &
      & units='W.m-2',                 &
      & xtype='float')

    call HistoryAutoAddVariable(                      &
      & varname='SfcXMomFlux',                        &
      & dims=(/'x','y','t'/),                         &
      & longname='surface x-component momentum flux', &
      & units='kg.m-2.s-1',                           &
      & xtype='float')

    call HistoryAutoAddVariable(                      &
      & varname='SfcYMomFlux',                        &
      & dims=(/'x','y','t'/),                         &
      & longname='surface y-component momentum flux', &
      & units='kg.m-2.s-1',                           &
      & xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable(                               &
        & varname=trim(SpcWetSymbol(l))//'_SfcMassFlux',         &
        & dims=(/'x','y','t'/),                              &
        & longname=trim(SpcWetSymbol(l))//' surface mass flux',  &
        & units='kg.m-2.s-1',                                    &
        & xtype='float')
    end do

  end subroutine Surfaceflux_Bulk_init
!!!------------------------------------------------------------------------!!!
  !---------------------------------------------------------------
  ! 地表面温度 my_sfcTemp を引数に追加 (mkuriki)
  ! subroutine Surfaceflux_Bulk_forcing( &
  !   &   pyz_VelX, xqz_VelY, xyz_PTemp, xyz_Exner, xyzf_QMix, &
  !   &   pyz_DVelXDt, xqz_DVelYDt, xyz_DPTempDt, xyz_DExnerDt, xyzf_DQMixDt &
  !   & )
  subroutine Surfaceflux_Bulk_forcing( &
    &   pyz_VelX, xqz_VelY, xyz_PTemp, xyz_Exner, xyzf_QMix, &
    &   pyz_DVelXDt, xqz_DVelYDt, xyz_DPTempDt, xyz_DExnerDt, xyzf_DQMixDt, &
    &   my_sfcTemp &
    & )
    ! 
    ! 下部境界からのフラックスによる温度の変化率を,
    ! バルク方法に基づいて計算する.
    !

    ! 暗黙の型宣言禁止
    ! Implicit none
    !
    implicit none

    ! 変数
    ! variables
    !
    real(DP), intent(in)   :: pyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速
                              ! X-component velocity
    real(DP), intent(in)   :: xqz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速
                              ! Y-component velocity
    real(DP), intent(in)   :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax)
                              ! 温位
                              ! Potential temperature
    real(DP), intent(in)   :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)
                              ! 圧力関数
                              ! Exner function
    real(DP), intent(in)   :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                              ! 混合比
                              ! Mixing ration
    real(DP), intent(inout):: pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax)
                              ! 風速時間変化率
                              ! X-component velocity tendency
    real(DP), intent(inout):: xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax)
                              ! 風速時間変化率
                              ! Y-component velocity tendency
    real(DP), intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax)
                              ! 温位時間変化率
                              ! Potential tempreture tendency
    real(DP), intent(inout):: xyz_DExnerDt(imin:imax,jmin:jmax,kmin:kmax)
                              ! 圧力関数時間変化率
                              ! Exner function tendency
    real(DP), intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                              ! 混合比時間変化率
                              ! Mixing ratio tendency

    ! 作業変数
    ! Work variables
    real(DP) :: xy_SurfRoughLength(imin:imax,jmin:jmax)
                              ! 祖度長さ
                              ! Roughness length
    real(DP) :: xy_SurfTemp(imin:imax,jmin:jmax)
                              ! 地表面温度
                              ! surface temperature
    real(DP) :: xy_SurfHumidCoef(imin:imax,jmin:jmax)
                              ! 
                              ! 
    real(DP) :: xy_SurfHeight(imin:imax,jmin:jmax)
                              ! 
                              ! 
    real(DP) :: xy_SurfVelTransCoef(imin:imax,jmin:jmax)
                              ! バルク係数(運動量)
                              ! Bulk coefficient for momentum
    real(DP) :: xy_SurfTempTransCoef(imin:imax,jmin:jmax)
                              ! バルク係数(熱)
                              ! Bulk coefficient for heat
    real(DP) :: xy_SurfQMixTransCoef(imin:imax,jmin:jmax)
                              ! バルク係数(混合比)
                              ! Bulk coefficient for mixing ratio
    real(DP) :: xyr_MomFluxX(imin:imax,jmin:jmax,kmin:kmax)
                              ! x 方向運動量フラックス
                              ! momentum flux in x
    real(DP) :: pyr_MomFluxX(imin:imax,jmin:jmax,kmin:kmax)
                              ! x 方向運動量フラックス
                              ! momentum flux in x
    real(DP) :: xyr_MomFluxY(imin:imax,jmin:jmax,kmin:kmax)
                              ! y 方向運動量フラックス
                              ! momentum flux in y
    real(DP) :: xqr_MomFluxY(imin:imax,jmin:jmax,kmin:kmax)
                              ! y 方向運動量フラックス
                              ! momentum flux in y
    real(DP) :: xyr_HeatFlux(imin:imax,jmin:jmax,kmin:kmax)
                              ! 熱フラックス
                              ! heat flux
    real(DP) :: xyrf_QMixFlux(imin:imax,jmin:jmax,kmin:kmax,ncmax)
                              ! 凝結成分混合比フラックス
                              ! Mixing ratio flux
    real(DP) :: xy_PTempFlux(imin:imax,jmin:jmax)
                              ! 温位フラックス
                              ! Potential temperature flux
    real(DP) :: py_VelXFlux(imin:imax,jmin:jmax)
                              ! 速度フラックス
                              ! velocity flux
    real(DP) :: xq_VelYFlux(imin:imax,jmin:jmax)
                              ! 速度フラックス
                              ! velocity flux
    real(DP) :: xyf_QMixFlux(imin:imax,jmin:jmax,ncmax)
                              ! 凝結成分混合比フラックス
                              ! Mixing ratio flux
    real(DP) :: xy_ExnerFlux(imin:imax,jmin:jmax)
                              !
                              !
    real(DP) :: xyz_VelX(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速 (xyz 格子)
                              ! X-component velocity (xyz grid)
    real(DP) :: xyz_VelY(imin:imax,jmin:jmax,kmin:kmax)
                              ! 水平風速 (xyz 格子)
                              ! Y-component velocity (xyz grid)
    real(DP) :: xyz_PTempAll(imin:imax,jmin:jmax,kmin:kmax)
                              ! 温位(基本場 + 擾乱)
                              ! Total value of potential temperature
    real(DP) :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax)
                              ! 温度(基本場 + 擾乱)
                              ! Total value of potential temperature
    real(DP) :: xyr_TempAll(imin:imax,jmin:jmax,kmin:kmax)
                              ! 温度(基本場 + 擾乱)
                              ! Total value of potential temperature
    real(DP) :: xyz_VirTempAll(imin:imax,jmin:jmax,kmin:kmax)
                              ! 仮温度(基本場 + 擾乱)
                              ! Total value of virtual temperature
    real(DP) :: xyr_VirTempAll(imin:imax,jmin:jmax,kmin:kmax)
                              ! 仮温度(基本場 + 擾乱)
                              ! Total value of virtual temperature
    real(DP) :: xyz_ExnerAll (imin:imax,jmin:jmax,kmin:kmax)
                              ! 圧力関数(基本場 + 擾乱)
                              ! Total value of exner function
    real(DP) :: xyr_ExnerAll (imin:imax,jmin:jmax,kmin:kmax)
                              ! 圧力関数(基本場 + 擾乱)
                              ! Total value of exner function
    real(DP) :: xyz_PressAll (imin:imax,jmin:jmax,kmin:kmax)
                              ! 圧力(基本場 + 擾乱)
                              ! Total value of pressure
    real(DP) :: xyr_PressAll (imin:imax,jmin:jmax,kmin:kmax)
                              ! 圧力(基本場 + 擾乱)
                              ! Total value of pressure
    real(DP) :: xyz_DensAll (imin:imax,jmin:jmax,kmin:kmax)
                              ! 圧力(基本場 + 擾乱)
                              ! Total value of pressure
    real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                              ! 混合比(基本場 + 擾乱)
                              ! Total value of mixing ratios
    real(DP) :: xyzf_QMixPerMolWt(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                              ! 混合比(基本場 + 擾乱)/分子量
                              ! mixing ratios per molecular weight

    real(DP) :: xy_DPTempDtBulk(imin:imax,jmin:jmax)
                              ! 時間変化(温位)
                              ! potential temperature tendency by surface flux
    real(DP) :: xy_DExnerDtBulk(imin:imax,jmin:jmax)
                              ! 時間変化(圧力関数)
                              ! Exner function tendency by surface flux
    real(DP) :: xyf_DQMixDtBulk(imin:imax,jmin:jmax, ncmax)
                              ! 時間変化(混合比)
                              ! Mixing ratio tendency by surface flux
    real(DP) :: py_DVelXDtBulk (imin:imax,jmin:jmax)
                              ! 時間変化(U)
                              ! x-component velocity tendency by surface flux
    real(DP) :: xq_DVelYDtBulk (imin:imax,jmin:jmax)
                              ! 時間変化(V)
                              ! y-component velocity tendency by surface flux
    real(DP) :: xyz_DPTempDtBulk(imin:imax,jmin:jmax,kmin:kmax)
                              ! 時間変化(温位)
                              ! potential temperature tendency by surface flux
    real(DP) :: xyz_DExnerDtBulk(imin:imax,jmin:jmax,kmin:kmax)
                              ! 時間変化(圧力関数)
                              ! Exner function tendency by surface flux
    real(DP) :: xyzf_DQMixDtBulk(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                              ! 時間変化(混合比)
                              ! Mixing ratio tendency by surface flux
    real(DP) :: pyz_DVelXDtBulk (imin:imax,jmin:jmax,kmin:kmax)
                              ! 時間変化(U)
                              ! x-component velocity tendency by surface flux
    real(DP) :: xqz_DVelYDtBulk (imin:imax,jmin:jmax,kmin:kmax)
                              ! 時間変化(V)
                              ! y-component velocity tendency by surface flux
    integer  :: kz            ! 配列添字
                              ! Arrzy index
    integer  :: s             ! 組成方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in dimension of constituen                         

    ! 初期化
    ! Initialization
    ! 
    kz = 1

!    xyr_MomFluxX = 0.0d0
!    xyr_MomFluxY = 0.0d0
!    xyr_HeatFlux = 0.0d0
!    xyrf_QMixFlux = 0.0d0
!    xy_SurfVelTransCoef = 0.0d0
!    xy_SurfTempTransCoef = 0.0d0
!    xy_SurfQMixTransCoef = 0.0d0

    ! 祖度長さの指定
    ! Specify surface length

    ! 表面値の指定
    ! Specify surface values
    xy_SurfHeight = 0.0d0
    ! ----------------------------------------
    ! 時刻ごとに変わる my_sfcTemp を地表面温度として代入 (mkuriki)
    ! xy_SurfTemp = TempSfc
    xy_SurfTemp = my_sfcTemp
    xy_SurfHumidCoef = 1.0d0
    xy_SurfRoughLength = SfcRoughLength

    ! 全量の計算
    ! Calculate total value of thermodynamic variables
    ! 
!    xyz_PTempAll  = xyz_PTemp + xyz_PTempBZ
!    xyz_ExnerAll  = xyz_Exner + xyz_ExnerBZ
!    xyzf_QMixAll  = xyzf_QMix + xyzf_QMixBZ

    xyz_PTempAll  = xyz_PTempBZ
    xyz_ExnerAll  = xyz_ExnerBZ
    xyzf_QMixAll  = xyzf_QMixBZ

    xyz_TempAll   = xyz_PTempAll * xyz_ExnerAll
    xyr_TempAll   = xyr_avr_xyz(xyz_TempAll)

    do s = 1, ncmax
      xyzf_QMixPerMolWt(:,:,:,s) = xyzf_QMixAll(:,:,:,IdxG(s)) / MolWtWet(IdxG(s))
    end do

!    xyz_VirTempAll = xyz_TempAll / &
!      & (                          &
!      &   (1.0d0 / ( 1.0d0 + MolWtDry * sum(xyzf_QMixPerMolWt, 4) ) ) /  & 
!      &   (1.0d0 + sum(xyzf_QMixAll, 4))     &
!      & )
!    xyr_VirTempAll = xyr_avr_xyz(xyz_VirTempAll)

    xyz_VirTempAll = xyz_TempAll
    xyr_VirTempAll = xyr_avr_xyz(xyz_TempAll)

    xyz_PressAll  = PressBasis * xyz_ExnerAll**(CpDry / GasRDry)
    xyr_PressAll  = xyr_avr_xyz(xyz_PressAll)
    xyr_ExnerAll  = xyr_avr_xyz(xyz_ExnerAll)

    xyz_DensAll   = xyz_PressAll / (GasRDry * xyz_VirTempAll)

    write(*,*) 'ExnerAll(0,0,0)=', xyr_ExnerAll(1,1,0) 
    write(*,*) 'ExnerAll(0,0,1)=', xyz_ExnerAll(1,1,1)
    write(*,*) 'TempAll(0,0,1)=', xyz_TempAll(1,1,1)
!    write(*,*) 'PTempAll(0,0,1)=', xyz_PTempAll(0,0,1)
    write(*,*) 'DensAll(0,0,1)=', xyz_DensAll(1,1,1)
    write(*,*) 'SurfTemp(0,0)=', xy_SurfTemp(1,1)
    

    ! Perturbation component of Exner function at the surface is assumed 
    ! to be same as that at the lowest layer. (YOT, 2011/09/03)
    ! 
    !ExnerBZSfc    = (PressSfc / PressBasis) ** (GasRDry / CpDry)
    !xy_PressSfc   = (PressBasis * xyz_ExnerAll(:,:,kz))**(CpDry / GasRDry)

    ! xyz 格子点の速度の計算
    ! Calculate velocities at xyz grid points
    !
    xyz_VelX = xyz_avr_pyz(pyz_VelX)
    xyz_VelY = xyz_avr_xqz(xqz_VelY)

    ! フラックスの計算
    ! Surface fluxes are calculated. 
    !
    ! dcpam ソースを call
    ! 
    xyrf_QMixFlux = 0.0_DP
    xyf_QMixFlux  = 0.0_DP

    call SurfaceFlux_Core( &
      & xyz_VelX   (1:nx,1:ny,1:nz), & ! (in)
      & xyz_VelY   (1:nx,1:ny,1:nz), & ! (in)
      & xyz_TempAll(1:nx,1:ny,1:nz), & ! (in)
      & xyr_TempAll(1:nx,1:ny,0:nz), & ! (in)
      & xyr_VirTempAll(1:nx,1:ny,0:nz), & ! (in)
      & xyzf_QMixAll(1:nx,1:ny,1:nz,1:ncmax), & ! (in)
      & xyr_PressAll(1:nx,1:ny,0:nz), & ! (in)
      & xy_SurfHeight(1:nx,1:ny),     & ! (in)
      & xyz_Z(1:nx,1:ny,1:nz),        & ! (in)
      & xyz_ExnerAll(1:nx,1:ny,1:nz), & ! (in)
      & xyr_ExnerAll(1:nx,1:ny,0:nz), & ! (in)
      & xy_SurfTemp(1:nx,1:ny),       & ! (in)
      & xy_SurfHumidCoef(1:nx,1:ny),  & ! (in)
      & xy_SurfRoughLength(1:nx,1:ny),& ! (in)
      & xyr_MomFluxX(1:nx,1:ny,0:nz), & ! (inout)
      & xyr_MomFluxY(1:nx,1:ny,0:nz), & ! (inout)
      & xyr_HeatFlux(1:nx,1:ny,0:nz), & ! (inout)
      & xyrf_QMixFlux(1:nx,1:ny,0:nz,1:ncmax), & ! (inout)
      & xy_SurfVelTransCoef(1:nx,1:ny), & ! (inout)
      & xy_SurfTempTransCoef(1:nx,1:ny),& ! (out)
      & xy_SurfQMixTransCoef(1:nx,1:ny) & ! (out)
      & )


    ! フラックスの変換
    ! convert surface flux
    !
    pyr_MomFluxX = pyr_avr_xyr(xyr_MomFluxX)
    xqr_MomFluxY = xqr_avr_xyr(xyr_MomFluxY)

    xy_PTempFlux(1:nx,1:ny) = xyr_HeatFlux(1:nx,1:ny,0)/ &
      &                         xyz_DensAll(1:nx,1:ny,kz) / CpDry
    py_VelXFlux(1:nx,1:ny)  = pyr_MomFluxX(1:nx,1:ny,0)/ &
      &                         xyz_DensAll(1:nx,1:ny,kz)
    xq_VelYFlux(1:nx,1:ny)  = xqr_MomFluxY(1:nx,1:ny,0)/ &
      &                         xyz_DensAll(1:nx,1:ny,kz)

    xyf_QMixFlux(1:nx,1:ny,:) = xyrf_QmixFlux(1:nx,1:ny,0,:)


    ! 圧力方程式の強制項の計算
    ! 
    xy_ExnerFlux = xy_DExnerDt_xy_xyf(xy_PTempFlux, xyf_QMixFlux, kz)


    ! 地表フラックスによる時間変化を計算
    ! Tendencies by surface fluxes (convergences of fluxes) are calculated.
    !
    xy_DPTempDtBulk = - ( 0.0d0 - xy_PTempFlux ) / z_dz(kz)
    !
    xy_DExnerDtBulk = - ( 0.0d0 - xy_ExnerFlux ) / z_dz(kz)
    !
    xyf_DQMixDtBulk = - ( 0.0d0 - xyf_QMixFlux ) / z_dz(kz)
    !
    py_DVelXDtBulk = - ( 0.0d0 - py_VelXFlux ) / z_dz(kz)
    !
    xq_DVelYDtBulk = - ( 0.0d0 - xq_VelYFlux ) / z_dz(kz)

    
    ! 仮引数配列へ格納
    ! Add tendency by surface flux convergence
    !
    xyz_DPTempDt(:,:,kz) = xyz_DPTempDt(:,:,kz) + xy_DPTempDtBulk
    xyz_DExnerDt(:,:,kz) = xyz_DExnerDt(:,:,kz) + xy_DExnerDtBulk
    do s = 1, ncmax
      xyzf_DQMixDt(:,:,kz,s) = xyzf_DQMixDt(:,:,kz,s) + xyf_DQMixDtBulk(:,:,s)
    end do
    pyz_DVelXDt (:,:,kz) = pyz_DVelXDt (:,:,kz) + py_DVelXDtBulk
    xqz_DVelYDt (:,:,kz) = xqz_DVelYDt (:,:,kz) + xq_DVelYDtBulk

    ! 出力
    ! Output
    !
    xyz_DPTempDtBulk = 0.0d0
    xyz_DExnerDtBulk = 0.0d0
    xyzf_DQMixDtBulk = 0.0d0
    pyz_DVelXDtBulk  = 0.0d0
    xqz_DVelYDtBulk  = 0.0d0

    xyz_DPTempDtBulk(:,:,kz) = xy_DPTempDtBulk
    xyz_DExnerDtBulk(:,:,kz) = xy_DExnerDtBulk
    do s = 1, ncmax
      xyzf_DQMixDtBulk(:,:,kz,s) = xyf_DQMixDtBulk(:,:,s)
    end do
    pyz_DVelXDtBulk (:,:,kz) = py_DVelXDtBulk
    xqz_DVelYDtBulk (:,:,kz) = xq_DVelYDtBulk
    !
    call HistoryAutoPut(TimeN, 'PTempSfc', xyz_DPTempDtBulk(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'ExnerSfc', xyz_DExnerDtBulk(1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelXSfc',  pyz_DVelXDtBulk (1:nx,1:ny,1:nz))
    call HistoryAutoPut(TimeN, 'VelYSfc',  xqz_DVelYDtBulk (1:nx,1:ny,1:nz))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_Sfc', &
        & xyzf_DQMixDtBulk(1:nx,1:ny,1:nz,s))
    end do

    call HistoryAutoPut(TimeN, 'PTempSfcFlux', xy_PTempFlux(1:nx,1:ny))
    call HistoryAutoPut(TimeN, 'ExnerSfcFlux', xy_ExnerFlux(1:nx,1:ny))
    call HistoryAutoPut(TimeN, 'VelXSfcFlux',  py_VelXFlux (1:nx,1:ny))
    call HistoryAutoPut(TimeN, 'VelYSfcFlux',  xq_VelYFlux (1:nx,1:ny))
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_SfcFlux', &
        & xyf_QMixFlux(1:nx,1:ny,s))
    end do

    call HistoryAutoPut(TimeN, 'SfcHeatFlux', xyr_HeatFlux(1:nx,1:ny,0) )
    call HistoryAutoPut(TimeN, 'SfcXMomFlux', xyr_MomFluxX(1:nx,1:ny,0) )
    call HistoryAutoPut(TimeN, 'SfcYMomFlux', xyr_MomFluxY(1:nx,1:ny,0) )
    do s = 1, ncmax
      call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_SfcMassFlux', &
        & xyz_DensAll(1:nx,1:ny,1) * xyf_QMixFlux(1:nx,1:ny,s))
    end do

  end subroutine Surfaceflux_Bulk_forcing
!!!------------------------------------------------------------------------!!!  
end module Surfaceflux_bulk
