subroutine PhySurfFluxBulkCoeff( phy_surfflx, xy_SurfBulkRiNum, xy_SurfVelRoughLength, xy_SurfTempRoughLength, xy_SurfGeoPot, xy_SurfVelBulkCoeff, xy_SurfTempBulkCoeff, xy_SurfQvapBulkCoeff, err )
    !
    ! バルク係数を算出します.
    !
    ! なお, 与えられた *phy_surfflx* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Bulk coefficients are calculated.
    !
    ! If *phy_surfflx* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYSURFFLX), intent(inout):: phy_surfflx
    real(DP), intent(in):: xy_SurfBulkRiNum (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $ number

!!$    real(DP), intent(in):: xy_SurfVelAbs (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
!!$                              ! 風速絶対値. 
!!$                              ! Absolute velocity
    real(DP), intent(in):: xy_SurfVelRoughLength (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 地表粗度長 (運動量). 
                              ! Surface rough length (momentum)
    real(DP), intent(in):: xy_SurfTempRoughLength (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 地表粗度長 (熱). 
                              ! Surface rough length (thermal)
    real(DP), intent(in):: xy_SurfGeoPot (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! $ \phi $ . ジオポテンシャル (地表面). 
                              ! Geo-potential (on surface)
    real(DP), intent(out):: xy_SurfVelBulkCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! バルク係数：運動量. 
                              ! Bulk coefficient: temperature
    real(DP), intent(out):: xy_SurfTempBulkCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! バルク係数：温度. 
                              ! Bulk coefficient: temperature
    real(DP), intent(out):: xy_SurfQvapBulkCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! バルク係数：比湿. 
                              ! Bulk coefficient: specific humidity
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: imax ! 経度格子点数. 
                   ! Number of grid points in longitude
    integer:: jmax ! 緯度格子点数. 
                   ! Number of grid points in latitude
    real(DP):: FKarm     ! $ k $ .                 カルマン定数.   Karman constant

    logical:: Neutral = .false.
                              ! 中立であるか否か. 
                              ! Whether it was neutral or not?
    real(DP), parameter:: ConstBulkCoeff  = -1.0_DP
                              ! バルク係数一定値. 
                              ! Steady value of bulk coefficient
    real(DP), parameter:: VelBulkCoeffMin  = 0.0_DP
                              ! $ u $ バルク係数最小値. 
                              ! Minimum value of $ u $ bulk coefficient
    real(DP), parameter:: TempBulkCoeffMin = 0.0_DP
                              ! $ T $ バルク係数最小値. 
                              ! Minimum value of $ T $ bulk coefficient
    real(DP), parameter:: QvapBulkCoeffMin = 0.0_DP
                              ! $ q $ バルク係数最小値. 
                              ! Minimum value of $ q $ bulk coefficient
    real(DP), parameter:: VelBulkCoeffMax  = 1.0_DP
                              ! $ u $ バルク係数最大値. 
                              ! Maximum value of $ u $ bulk coefficient
    real(DP), parameter:: TempBulkCoeffMax = 1.0_DP
                              ! $ T $ バルク係数最大値. 
                              ! Maximum value of $ T $ bulk coefficient
    real(DP), parameter:: QvapBulkCoeffMax = 1.0_DP
                              ! $ q $ バルク係数最大値. 
                              ! Maximum value of $ q $ bulk coefficient

    integer:: i, j            ! DO ループ用作業変数
                              ! Work variables for DO loop
    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhySurfFluxBulkCoeff'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_surfflx % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYSURFFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  *phy_surfflx* に格納されている設定値の取り出し
    !  Fetch setting values stored in *phy_surfflx*
    !-----------------------------------------------------------------
    imax  = phy_surfflx % imax 
    jmax  = phy_surfflx % jmax 

    FKarm  = phy_surfflx % FKarm 

    !-----------------------------------------------------------------
    !  中立バルク係数の計算
    !  Calculate neutral bulk coefficients
    !-----------------------------------------------------------------
    if ( ConstBulkCoeff < 0.0_DP ) then
      
      xy_SurfVelBulkCoeff  = ( FKarm / log ( xy_SurfGeoPot / xy_SurfVelRoughLength ) )**2

      xy_SurfTempBulkCoeff = ( FKarm / log ( xy_SurfGeoPot / xy_SurfTempRoughLength ) )**2

      xy_SurfQvapBulkCoeff = xy_SurfTempBulkCoeff
      
    else
      xy_SurfVelBulkCoeff  = ConstBulkCoeff
      xy_SurfTempBulkCoeff = ConstBulkCoeff
      xy_SurfQvapBulkCoeff = ConstBulkCoeff
    end if
    
    !-----------------------------------------------------------------
    !  非中立バルク係数の計算
    !  Calculate non-neutral bulk coefficients
    !-----------------------------------------------------------------
    if ( .not. Neutral ) then
      
      do i = 0, imax-1
        do j = 0, jmax-1

          if ( xy_SurfBulkRiNum(i,j) > 0.0_DP ) then 

            xy_SurfVelBulkCoeff(i,j) = xy_SurfVelBulkCoeff(i,j) / (   1.0_DP + 10.0_DP * xy_SurfBulkRiNum(i,j) / sqrt( 1.0_DP + 5.0_DP * xy_SurfBulkRiNum(i,j) ) )

            xy_SurfTempBulkCoeff(i,j) = xy_SurfTempBulkCoeff(i,j) / (   1.0_DP + 15.0_DP * xy_SurfBulkRiNum(i,j) / sqrt( 1.0_DP + 5.0_DP * xy_SurfBulkRiNum(i,j) ) )

            xy_SurfQvapBulkCoeff(i,j) = xy_SurfTempBulkCoeff(i,j)

          else

            xy_SurfVelBulkCoeff(i,j) = xy_SurfVelBulkCoeff(i,j) * (   1.0_DP - 10.0_DP * xy_SurfBulkRiNum(i,j) / (   1.0_DP + 75.0_DP * xy_SurfVelBulkCoeff(i,j) * sqrt( - xy_SurfGeoPot(i,j) / xy_SurfVelRoughLength(i,j) * xy_SurfBulkRiNum(i,j) ) ) )
            
            xy_SurfTempBulkCoeff(i,j) = xy_SurfTempBulkCoeff(i,j) * (   1.0_DP - 15.0_DP * xy_SurfBulkRiNum(i,j) / (   1.0_DP + 75.0_DP * xy_SurfTempBulkCoeff(i,j) * sqrt( - xy_SurfGeoPot(i,j) / xy_SurfTempRoughLength(i,j) * xy_SurfBulkRiNum(i,j) ) ) )

            xy_SurfQvapBulkCoeff(i,j) = xy_SurfTempBulkCoeff(i,j)

          end if
        end do
      end do
      
    end if
    
    !-----------------------------------------------------------------
    !  最大/最小 判定
    !  Measure maximum/minimum
    !-----------------------------------------------------------------
    do i = 0, imax-1
      do j = 0, jmax-1

        xy_SurfVelBulkCoeff(i,j)  = max( min( xy_SurfVelBulkCoeff(i,j), VelBulkCoeffMax ), VelBulkCoeffMin )

        xy_SurfTempBulkCoeff(i,j) = max( min( xy_SurfTempBulkCoeff(i,j), TempBulkCoeffMax ), TempBulkCoeffMin )

        xy_SurfQvapBulkCoeff(i,j) = max( min( xy_SurfQvapBulkCoeff(i,j), QvapBulkCoeffMax ), QvapBulkCoeffMin )

      end do
    end do
    
    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhySurfFluxBulkCoeff

  subroutine PhySurfFluxClose( phy_surfflx, err )
    !
    ! PHYSURFFLX 型の変数の終了処理を行います.
    ! なお, 与えられた *phy_surfflx* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Deconstructor of "PHYSURFFLX".
    ! Note that if *phy_surfflx* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYSURFFLX), intent(inout):: phy_surfflx
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhySurfFluxClose'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_surfflx % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYSURFFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  "PHYSURFFLX" の設定の消去
    !  Clear the settings for "PHYSURFFLX"
    !-----------------------------------------------------------------
    deallocate( phy_surfflx % xy_SurfTemp )
    deallocate( phy_surfflx % xy_SurfHumidCoeff )
    deallocate( phy_surfflx % xy_SurfRoughLength )
    deallocate( phy_surfflx % xy_SurfCondition )

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
    phy_surfflx % initialized = .false.
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhySurfFluxClose

  subroutine PhySurfFluxCreate( phy_surfflx, imax, jmax, kmax, RAir, Cp, Grav, EL, ES0, RVap, EpsV, FKarm, xy_SurfTemp, xy_SurfHumidCoeff, xy_SurfRoughLength, xy_SurfCondition, nmlfile, err )
    !
    ! PHYSURFFLX 型の変数の初期設定を行います.
    ! 他のサブルーチンを使用する前に必ずこのサブルーチンによって
    ! PHYSURFFLX 型の変数を初期設定してください.
    !
    ! なお, 与えられた *phy_surfflx* が既に初期設定されている場合,
    ! プログラムはエラーを発生させます.
    !
    ! NAMELIST を利用する場合には引数 *nmlfile* に NAMELIST ファイル名
    ! を与えてください. NAMELIST 変数群の詳細に関しては 
    ! NAMELIST#phy_surface_flux_nml を参照してください. 
    !
    ! Constructor of "PHYSURFFLX".
    ! Initialize *phy_surfflx* by this subroutine, 
    ! before other procedures are used, 
    !
    ! Note that if *phy_surfflx* is already initialized 
    ! by this procedure, error is occurred.
    !
    ! In order to use NAMELIST, specify a NAMELIST filename to 
    ! argument *nmlfile*. See "NAMELIST#phy_surface_flux_nml"
    ! for details about a NAMELIST group.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_present, only: present_and_not_empty, present_and_true
    use dc_message, only: MessageNotify
    use dc_error, only: StoreError, DC_NOERR, DC_EALREADYINIT, DC_EARGLACK, DC_ENEGATIVE, DC_ENOFILEREAD
    implicit none
    type(PHYSURFFLX), intent(inout):: phy_surfflx
    integer, intent(in):: imax ! 経度格子点数. 
                   ! Number of grid points in longitude
    integer, intent(in):: jmax ! 緯度格子点数. 
                   ! Number of grid points in latitude
    integer, intent(in):: kmax ! 鉛直層数. 
                   ! Number of vertical level
    real(DP), intent(in):: Grav      ! $ g $ .      重力加速度.     Gravitational acceleration
    real(DP), intent(in):: Cp        ! $ C_p $ .    大気定圧比熱.   Specific heat of air at constant pressure
    real(DP), intent(in):: RAir      ! $ R $ .      大気気体定数.   Gas constant of air
    real(DP), intent(in):: EL        ! $ L $ .      水の凝結の潜熱. Latent heat of condensation of water vapor
    real(DP), intent(in):: RVap      ! $ R_v $ .    水蒸気気体定数. Gas constant of water vapor
    real(DP), intent(in):: EpsV      ! $ \epsilon_v $ .        水蒸気分子量比. Molecular weight of water vapor
    real(DP), intent(in):: ES0       ! $ e^{*} $ (273K) .      0 ℃での飽和蒸気圧. Saturated vapor pressure at 0 degrees C
    real(DP), intent(in):: FKarm     ! $ k $ .                 カルマン定数.   Karman constant
    real(DP), intent(in):: xy_SurfTemp (0:imax-1, 0:jmax-1)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP), intent(in):: xy_SurfHumidCoeff (0:imax-1, 0:jmax-1)
                              ! 地表湿潤度. 
                              ! Surface humidity coefficient
    real(DP), intent(in):: xy_SurfRoughLength (0:imax-1, 0:jmax-1)
                              ! 地表粗度長. 
                              ! Surface rough length
    integer, intent(in):: xy_SurfCondition (0:imax-1, 0:jmax-1)
                              ! 地表状態. 
                              ! Surface condition
    character(*), intent(in), optional:: nmlfile
                              ! NAMELIST ファイルの名称. 
                              ! この引数に空文字以外を与えた場合, 
                              ! 指定されたファイルから 
                              ! NAMELIST 変数群を読み込みます. 
                              ! ファイルを読み込めない場合にはエラーを
                              ! 生じます.
                              !
                              ! NAMELIST 変数群の詳細に関しては 
                              ! NAMELIST#phy_surface_flux_nml 
                              ! を参照してください. 
                              !
                              ! NAMELIST file name. 
                              ! If nonnull character is specified to
                              ! this argument, 
                              ! NAMELIST group name is loaded from the 
                              ! file. 
                              ! If the file can not be read, 
                              ! an error occurs.
                              ! 
                              ! See "NAMELIST#phy_surface_flux_nml" 
                              ! for details about a NAMELIST group.
                              ! 
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhySurfFluxCreate'
  continue
    call BeginSub( subname, version )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( phy_surfflx % initialized ) then
      stat = DC_EALREADYINIT
      cause_c = 'PHYSURFFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  引数の正当性のチェック
    !  Validate arguments
    !-----------------------------------------------------------------
    if (imax < 1) then
      stat = DC_ENEGATIVE
      cause_c = 'imax'
      goto 999
    end if
    if (jmax < 1) then
      stat = DC_ENEGATIVE
      cause_c = 'jmax'
      goto 999
    end if
    if (kmax < 1) then
      stat = DC_ENEGATIVE
      cause_c = 'kmax'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  "PHYSURFFLX" の設定
    !  Configure the settings for "PHYSURFFLX"
    !-----------------------------------------------------------------
    phy_surfflx % imax  = imax 
    phy_surfflx % jmax  = jmax 
    phy_surfflx % kmax  = kmax 
    phy_surfflx % RAir  = RAir 
    phy_surfflx % Cp    = Cp   
    phy_surfflx % Grav  = Grav 
    phy_surfflx % EL    = EL   
    phy_surfflx % ES0   = ES0  
    phy_surfflx % RVap  = RVap 
    phy_surfflx % EpsV  = EpsV 
    phy_surfflx % FKarm = FKarm

    allocate( phy_surfflx % xy_SurfTemp (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1) )
    allocate( phy_surfflx % xy_SurfHumidCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1) )
    allocate( phy_surfflx % xy_SurfRoughLength (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1) )
    allocate( phy_surfflx % xy_SurfCondition (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1) )

    phy_surfflx % xy_SurfTemp        = xy_SurfTemp       
    phy_surfflx % xy_SurfHumidCoeff  = xy_SurfHumidCoeff 
    phy_surfflx % xy_SurfRoughLength = xy_SurfRoughLength
    phy_surfflx % xy_SurfCondition   = xy_SurfCondition  

    !-------------------------
    !  デフォルト値
    !  Default values
!!$    phy_surfflx % param_r = 0.0_DP
!!$    phy_surfflx % param_c = 'hogehoge'

    !-------------------------
    !  オプショナル引数からの値
    !  Values from optional arguments
!!$    phy_surfflx % param_i = param_i
!!$    if ( present(param_r) )  phy_surfflx % param_r = param_r
!!$    if ( present(param_c) )  phy_surfflx % param_c = param_c

    !-------------------------
    !  NAMELIST からの値
    !  Values from NAMELIST

!!$    if ( present_and_not_empty(nmlfile) ) then
!!$      call MessageNotify( 'M', subname, &
!!$        & 'Loading NAMELIST file "%c" ...', &
!!$        & c1 = trim(nmlfile) )
!!$      call NmlRead ( nmlfile = nmlfile, &      ! (in)
!!$        & param_i = phy_surfflx % param_i, &   ! (inout)
!!$        & param_r = phy_surfflx % param_r, &   ! (inout)
!!$        & param_c_ = phy_surfflx % param_c, &  ! (inout)
!!$        & err = err )                          ! (out)
!!$      if ( present_and_true(err) ) then
!!$        call MessageNotify( 'W', subname, &
!!$          & '"%c" can not be read.', &
!!$          & c1 = trim(nmlfile) )
!!$        stat = DC_ENOFILEREAD
!!$        cause_c = nmlfile
!!$        goto 999
!!$      end if
!!$    end if

    !-----------------------------------------------------------------
    !  設定値の正当性のチェック
    !  Validate setting values
    !-----------------------------------------------------------------
!!$    if ( phy_surfflx % param_i < 0 ) then
!!$      stat = DC_ENEGATIVE
!!$      cause_c = 'param_i'
!!$      goto 999
!!$    end if


    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
    phy_surfflx % initialized = .true.
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhySurfFluxCreate

  logical function PhySurfFluxInitialized( phy_surfflx ) result(result)
    !
    ! *phy_surfflx* が初期設定されている場合には .true. が,
    ! 初期設定されていない場合には .false. が返ります.
    !
    ! If *phy_surfflx* is initialized, .true. is returned.
    ! If *phy_surfflx* is not initialized, .false. is returned.
    !
    implicit none
    type(PHYSURFFLX), intent(in):: phy_surfflx
  continue
    result = phy_surfflx % initialized
  end function PhySurfFluxInitialized

  subroutine PhySurfFluxNmlRead( nmlfile, err )
    !
    ! NAMELIST ファイル *nmlfile* から値を入力するための
    ! 内部サブルーチンです. Create 内で呼び出されることを
    ! 想定しています.
    !
    ! 値が NAMELIST ファイル内で指定されていない場合には,
    ! 入力された値がそのまま返ります.
    !
    ! なお, *nmlfile* に空文字が与えられた場合, または
    ! 与えられた *nmlfile* を読み込むことができない場合, 
    ! プログラムはエラーを発生させます.
    !
    ! This is an internal subroutine to input values from 
    ! NAMELIST file *nmlfile*. This subroutine is expected to be
    ! called by "Create".
    !
    ! A value not specified in NAMELIST file is returned
    ! without change.
    !
    ! If *nmlfile* is empty, or *nmlfile* can not be read, 
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_iounit, only: FileOpen
    use dc_message, only: MessageNotify
    use dc_present, only: present_and_true
    use dc_error, only: StoreError, DC_NOERR, DC_ENOFILEREAD
    implicit none
    character(*), intent(in):: nmlfile
                              ! NAMELIST ファイルの名称. 
                              ! NAMELIST file name
!!$    integer, intent(inout):: param_i
!!$    real(DP), intent(inout):: param_r
!!$    character(*), intent(inout):: param_c_
!!$    character(TOKEN):: param_c
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

!!$    namelist /phy_surface_flux_nml/ &
!!$      & param_i, param_r, param_c
                              ! phy_surface_flux モジュール用
                              ! NAMELIST 変数群名.
                              !
                              ! phy_surface_flux#Create を使用する際に, 
                              ! オプショナル引数 *nmlfile* へ NAMELIST 
                              ! ファイル名を指定することで, そのファイルから
                              ! この NAMELIST 変数群を読み込みます.
                              !
                              ! NAMELIST group name for 
                              ! "phy_surface_flux" module.
                              ! 
                              ! If a NAMELIST filename is specified to 
                              ! an optional argument *nmlfile* 
                              ! when "phy_surface_flux#Create" is used, 
                              ! this NAMELIST group is loaded from 
                              ! the file.

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号. 
                              ! Unit number for NAMELIST file open
!!$    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT. 
!!$                              ! IOSTAT of NAMELIST read
    character(*), parameter:: subname = 'PhySurfFluxNmlRead'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''



    !-----------------------------------------------------------------
    !  文字型引数を NAMELIST 変数群へ代入
    !  Substitute character arguments to NAMELIST group
    !-----------------------------------------------------------------
!!$    param_c = param_c_

    !----------------------------------------------------------------
    !  NAMELIST ファイルのオープン
    !  Open NAMELIST file
    !----------------------------------------------------------------
    call FileOpen( unit = unit_nml, file = nmlfile, mode = 'r', err = err )                   ! (out)
    if ( present_and_true(err) ) then
      stat = DC_ENOFILEREAD
      cause_c = nmlfile
      goto 999
    end if


    !-----------------------------------------------------------------
    !  NAMELIST 変数群の取得
    !  Get NAMELIST group
    !-----------------------------------------------------------------
!!$    read( unit = unit_nml, & ! (in)
!!$      & nml = phy_surface_flux_nml, iostat = iostat_nml ) ! (out)
!!$    if ( iostat_nml == 0 ) then
!!$      call MessageNotify( 'M', subname, &
!!$        & 'NAMELIST group "%c" is loaded from "%c".', &
!!$        & c1 = 'phy_surface_flux_nml', c2 = trim(nmlfile) )
!!$      write(STDOUT, nml = phy_surface_flux_nml)
!!$    else
!!$      call MessageNotify( 'W', subname, &
!!$        & 'NAMELIST group "%c" is not found in "%c" (iostat=%d).', &
!!$        & c1 = 'phy_surface_flux_nml', c2 = trim(nmlfile), &
!!$        & i = (/iostat_nml/) )
!!$    end if

    close( unit_nml )

    !-----------------------------------------------------------------
    !  NAMELIST 変数群を文字型引数へ代入
    !  Substitute NAMELIST group to character arguments
    !-----------------------------------------------------------------
!!$    param_c_ = param_c

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhySurfFluxNmlRead

  subroutine PhySurfFluxPutLine( phy_surfflx, unit, indent, err )
    !
    ! 引数 *phy_surfflx* に設定されている情報を印字します.
    ! デフォルトではメッセージは標準出力に出力されます. 
    ! *unit* に装置番号を指定することで, 出力先を変更することが可能です.
    !
    ! Print information of *phy_surfflx*.
    ! By default messages are output to standard output.
    ! Unit number for output can be changed by *unit* argument.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYSURFFLX), intent(in):: phy_surfflx
    integer, intent(in), optional:: unit
                              ! 出力先の装置番号.
                              ! デフォルトの出力先は標準出力.
                              !
                              ! Unit number for output.
                              ! Default value is standard output.
    character(*), intent(in), optional:: indent
                              ! 表示されるメッセージの字下げ.
                              !
                              ! Indent of displayed messages.
    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: stat
    character(STRING):: cause_c
    integer:: out_unit
    integer:: indent_len
    character(STRING):: indent_str
    character(*), parameter:: subname = 'PhySurfFluxPutLine'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( present(unit) ) then
      out_unit = unit
    else
      out_unit = STDOUT
    end if

    indent_len = 0
    indent_str = ''
    if ( present(indent) ) then
      if ( len(indent) /= 0 ) then
        indent_len = len(indent)
        indent_str(1:indent_len) = indent
      end if
    end if


    !-----------------------------------------------------------------
    !  "PHYSURFFLX" の設定の印字
    !  Print the settings for "PHYSURFFLX"
    !-----------------------------------------------------------------
    if ( phy_surfflx % initialized ) then
      call Printf( out_unit, indent_str(1:indent_len) // '#<PHYSURFFLX:: @initialized=%y', l = (/phy_surfflx % initialized/) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @imax=%d @jmax=%d @kmax=%d', i = (/phy_surfflx % imax, phy_surfflx % jmax, phy_surfflx % kmax/) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @Grav=%f @Cp=%f @RAir=%f', d = (/ phy_surfflx % Grav, phy_surfflx % Cp, phy_surfflx % RAir /) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @EL=%f @RVap=%f @EpsV=%f', d = (/ phy_surfflx % EL, phy_surfflx % RVap, phy_surfflx % EpsV /) )

      call Printf( out_unit, indent_str(1:indent_len) // ' @ES0=%f @FKarm=%f', d = (/ phy_surfflx % ES0, phy_surfflx % FKarm /) )

      call PutLine( phy_surfflx % xy_SurfTemp, unit = out_unit, lbounds = lbound(phy_surfflx % xy_SurfTemp), ubounds = ubound(phy_surfflx % xy_SurfTemp), indent = indent_str(1:indent_len) // ' @xy_SurfTemp=' )

      call PutLine( phy_surfflx % xy_SurfHumidCoeff, unit = out_unit, lbounds = lbound(phy_surfflx % xy_SurfHumidCoeff), ubounds = ubound(phy_surfflx % xy_SurfHumidCoeff), indent = indent_str(1:indent_len) // ' @xy_SurfHumidCoeff=' )

      call PutLine( phy_surfflx % xy_SurfRoughLength, unit = out_unit, lbounds = lbound(phy_surfflx % xy_SurfRoughLength), ubounds = ubound(phy_surfflx % xy_SurfRoughLength), indent = indent_str(1:indent_len) // ' @xy_SurfRoughLength=' )

      call PutLine( phy_surfflx % xy_SurfCondition, unit = out_unit, lbounds = lbound(phy_surfflx % xy_SurfCondition), ubounds = ubound(phy_surfflx % xy_SurfCondition), indent = indent_str(1:indent_len) // ' @xy_SurfCondition=' )

      call Printf( out_unit, indent_str(1:indent_len) // '>' )
    else
      call Printf( out_unit, indent_str(1:indent_len) // '#<PHYSURFFLX:: @initialized=%y>', l = (/phy_surfflx % initialized/) )
    end if

    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhySurfFluxPutLine

  subroutine PhySurfFluxSurfaceFlux( phy_surfflx, xyz_U, xyz_V, xyz_Temp, xyr_Temp, xyz_QVap, xyz_Press, xyr_Press, xyz_GeoPot, xyr_UFlux, xyr_VFlux, xyr_TempFlux, xyr_QVapFlux, xy_SurfUVMatrix, xyaa_SurfTempMatrix, xyaa_SurfQVapMatrix, err )
    !
    ! 地表面フラックスを計算します.
    !
    ! なお, 与えられた *phy_surfflx* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Surface fluxes are calculated.
    !
    ! If *phy_surfflx* is not initialized by "Create" yet,
    ! error is occurred.
    !
    use dc_trace, only: BeginSub, EndSub
    use dc_string, only: PutLine, Printf
    use dc_types, only: DP, STRING, TOKEN, STDOUT
    use dc_error, only: StoreError, DC_NOERR, DC_ENOTINIT
    implicit none
    type(PHYSURFFLX), intent(inout):: phy_surfflx
    real(DP), intent(in):: xyz_U (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax-1)
                              ! $ U $ . 東西風速. Zonal wind
    real(DP), intent(in):: xyz_V (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax-1)
                              ! $ V $ . 南北風速. Meridional wind

    real(DP), intent(in):: xyz_Temp (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax-1)
                              ! $ T $ . 温度 (整数レベル). 
                              ! Temperature (full level)
    real(DP), intent(in):: xyr_Temp (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax)
                              ! $ T $ . 温度 (半整数レベル). 
                              ! Temperature (half level)
    real(DP), intent(in):: xyz_QVap (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax-1)
                              ! $ q $ .     比湿. Specific humidity
    real(DP), intent(in):: xyz_Press (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax-1)
                              ! $ P_s $ . 地表面気圧 (整数レベル). 
                              ! Surface pressure (full level)
    real(DP), intent(in):: xyr_Press (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax)
                              ! $ P_s $ . 地表面気圧 (半整数レベル). 
                              ! Surface pressure (half level)
    real(DP), intent(in):: xyz_GeoPot (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax-1)
                              ! $ \phi $ . ジオポテンシャル (整数レベル). 
                              ! Geo-potential (full level)
    real(DP), intent(inout):: xyr_UFlux (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax)
                              ! 東西風速フラックス. 
                              ! Zonal wind flux
    real(DP), intent(inout):: xyr_VFlux (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax)
                              ! 南北風速フラックス. 
                              ! Meridional wind flux
    real(DP), intent(inout):: xyr_TempFlux (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax)
                              ! 温度フラックス. 
                              ! Temperature flux
    real(DP), intent(inout):: xyr_QvapFlux (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:phy_surfflx%kmax)
                              ! 比湿フラックス. 
                              ! Specific humidity flux
    real(DP), intent(out):: xy_SurfUVMatrix (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              !  速度陰解行列: 地表. 
                              ! Implicit matrix about velocity: surface
    real(DP), intent(out):: xyaa_SurfTempMatrix (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:1, -1:1)

                              ! 温度陰解行列: 地表. 
                              ! Implicit matrix about temperature: surface
    real(DP), intent(out):: xyaa_SurfQvapMatrix (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1, 0:1, -1:1)
                              ! 比湿陰解行列: 地表. 
                              ! Implicit matrix about specific humidity: surface

    logical, intent(out), optional:: err
                              ! 例外処理用フラグ.
                              ! デフォルトでは, この手続き内でエラーが
                              ! 生じた場合, プログラムは強制終了します.
                              ! 引数 *err* が与えられる場合,
                              ! プログラムは強制終了せず, 代わりに
                              ! *err* に .true. が代入されます.
                              !
                              ! Exception handling flag. 
                              ! By default, when error occur in 
                              ! this procedure, the program aborts. 
                              ! If this *err* argument is given, 
                              ! .true. is substituted to *err* and 
                              ! the program does not abort. 

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: imax ! 経度格子点数. 
                   ! Number of grid points in longitude
    integer:: jmax ! 緯度格子点数. 
                   ! Number of grid points in latitude

    real(DP):: Grav      ! $ g $ .      重力加速度.     Gravitational acceleration
    real(DP):: Cp        ! $ C_p $ .    大気定圧比熱.   Specific heat of air at constant pressure
    real(DP):: RAir      ! $ R $ .      大気気体定数.   Gas constant of air
    real(DP):: EL        ! $ L $ .      水の凝結の潜熱. Latent heat of condensation of water vapor
    real(DP):: RVap      ! $ R_v $ .    水蒸気気体定数. Gas constant of water vapor
    real(DP):: EpsV      ! $ \epsilon_v $ .        水蒸気分子量比. Molecular weight of water vapor
    real(DP):: ES0       ! $ e^{*} $ (273K) .      0 ℃での飽和蒸気圧. Saturated vapor pressure at 0 degrees C

    real(DP):: xy_SurfTemp (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 地表面温度. 
                              ! Surface temperature
    real(DP):: xy_SurfHumidCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 地表湿潤度. 
                              ! Surface humidity coefficient
    real(DP):: xy_SurfRoughLength (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 地表粗度長. 
                              ! Surface rough length
    integer:: xy_SurfCondition (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 地表状態. 
                              ! Surface condition

    real(DP), parameter:: BasePotTemp   = 300.0_DP
                              ! 基本温位. 
                              ! Basic potential temperature
    real(DP), parameter:: VelMinForRi   = 0.01_DP
                              ! $ R_i $ 数用風最小値. 
                              ! Minimum value of velocity for $ R_i $ number
    real(DP), parameter:: VelMinForVel  = 0.01_DP
                              ! 運動量用風最小値. 
                              ! Minimum value of velocity for momentum
    real(DP), parameter:: VelMinForTemp = 0.01_DP
                              ! 熱用風最小値. 
                              ! Minimum value of velocity for thermal
    real(DP), parameter:: VelMinForQvap = 0.01_DP
                              ! 水蒸気用風最小値. 
                              ! Minimum value of velocity for vapor
    real(DP), parameter:: VelMaxForVel  = 1000.0_DP
                              ! 運動量用風最大値. 
                              ! Maximum value of velocity for momentum
    real(DP), parameter:: VelMaxForTemp = 1000.0_DP
                              ! 熱用風最大値. 
                              ! Maximum value of velocity for thermal
    real(DP), parameter:: VelMaxForQvap = 1000.0_DP
                              ! 水蒸気用風最大値. 
                              ! Maximum value of velocity for vapor

    real(DP):: xy_SurfBulkRiNum (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $ number
    real(DP):: xy_SurfTempTransCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP):: xy_SurfQvapTransCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 輸送係数：比湿. 
                              ! Transfer coefficient: specific humidity
    real(DP):: xy_SurfVelTransCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 輸送係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP):: xy_SurfTempBulkCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! バルク係数：温度. 
                              ! Bulk coefficient: temperature
    real(DP):: xy_SurfQvapBulkCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! バルク係数：比湿. 
                              ! Bulk coefficient: specific humidity
    real(DP):: xy_SurfVelBulkCoeff (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! バルク係数：運動量. 
                              ! Bulk coefficient: temperature
    real(DP):: xy_SurfExner (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! Exner 関数. 
                              ! Exner function
    real(DP):: xy_SurfVelAbs (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 風速絶対値. 
                              ! Absolute velocity
    real(DP):: xy_SurfSatQvap (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 地表飽和比湿. 
                              ! Saturated specific humidity on surface
    real(DP):: xy_SurfDSatQvapDTemp (0:phy_surfflx%imax-1, 0:phy_surfflx%jmax-1)
                              ! 地表飽和比湿変化. 
                              ! Saturated specific humidity tendency on surface

    integer:: i, j            ! DO ループ用作業変数
                              ! Work variables for DO loop
    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhySurfFluxSurfaceFlux'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( .not. phy_surfflx % initialized ) then
      stat = DC_ENOTINIT
      cause_c = 'PHYSURFFLX'
      goto 999
    end if

    !-----------------------------------------------------------------
    !  *phy_surfflx* に格納されている設定値の取り出し
    !  Fetch setting values stored in *phy_surfflx*
    !-----------------------------------------------------------------
    imax  = phy_surfflx % imax 
    jmax  = phy_surfflx % jmax 
    RAir  = phy_surfflx % RAir 
    Cp    = phy_surfflx % Cp   
    Grav  = phy_surfflx % Grav 
    EL    = phy_surfflx % EL   
    ES0   = phy_surfflx % ES0  
    RVap  = phy_surfflx % RVap 
    EpsV  = phy_surfflx % EpsV 

    xy_SurfTemp        = phy_surfflx % xy_SurfTemp       
    xy_SurfHumidCoeff  = phy_surfflx % xy_SurfHumidCoeff 
    xy_SurfRoughLength = phy_surfflx % xy_SurfRoughLength
    xy_SurfCondition   = phy_surfflx % xy_SurfCondition  


    !-----------------------------------------------------------------
    !  Exner 関数算出
    !  Calculate Exner functions
    !-----------------------------------------------------------------
    xy_SurfExner = ( xyz_Press(:,:,0) / xyr_Press(:,:,0)  )**( RAir / Cp )

    !-----------------------------------------------------------------
    !  バルク $ R_i $ 数算出
    !  Calculate bulk $ R_i $
    !-----------------------------------------------------------------
    do i = 0, imax-1
      do j = 0, jmax-1
        xy_SurfVelAbs(i,j) = sqrt ( xyz_U(i,j,0)**2 + xyz_V(i,j,0)**2 )

        xy_SurfBulkRiNum(i,j) = Grav / BasePotTemp * ( xyz_Temp(i,j,0) / xy_SurfExner(i,j) - xy_SurfTemp(i,j)  ) / max( xy_SurfVelAbs(i,j), VelMinForRi )**2 * xyz_GeoPot(i,j,0)
      end do
    end do
    
    !-----------------------------------------------------------------
    !  バルク係数算出
    !  Calculate bulk coefficients
    !-----------------------------------------------------------------
    call BulkCoeff( phy_surfflx = phy_surfflx, xy_SurfBulkRiNum       = xy_SurfBulkRiNum, xy_SurfVelRoughLength  = xy_SurfRoughLength, xy_SurfTempRoughLength = xy_SurfRoughLength, xy_SurfGeoPot          = xyz_GeoPot(:,:,0), xy_SurfVelBulkCoeff    = xy_SurfVelBulkCoeff, xy_SurfTempBulkCoeff   = xy_SurfTempBulkCoeff, xy_SurfQvapBulkCoeff   = xy_SurfQvapBulkCoeff )  ! (out)

    !-----------------------------------------------------------------
    !  輸送係数の計算
    !  Calculate transfer coefficient
    !-----------------------------------------------------------------
    do i = 0, imax-1
      do j = 0, jmax-1
        
        xy_SurfVelTransCoeff(i,j) = xy_SurfVelBulkCoeff(i,j) * xyr_Press(i,j,0) / ( RAir * xyr_Temp(i,j,0) ) * min( max( xy_SurfVelAbs(i,j), VelMinForVel ), VelMaxForVel )
        
        xy_SurfTempTransCoeff(i,j) = xy_SurfTempBulkCoeff(i,j) * xyr_Press(i,j,0) / ( RAir * xyr_Temp(i,j,0) ) * min( max( xy_SurfVelAbs(i,j), VelMinForTemp ), VelMaxForTemp )
        
        xy_SurfQvapTransCoeff(i,j) = xy_SurfQvapBulkCoeff(i,j) * xyr_Press(i,j,0) / ( RAir * xyr_Temp(i,j,0) ) * min( max( xy_SurfVelAbs(i,j), VelMinForQvap ), VelMaxForQvap )
        
      end do
    end do
    
    !-----------------------------------------------------------------
    !  飽和比湿の計算
    !  Calculate saturated specific humidity
    !-----------------------------------------------------------------
    xy_SurfSatQvap = EpsV * ES0 * exp( EL / RVap * ( 1.0_DP / 273.0_DP - 1.0_DP / xy_SurfTemp ) ) / xyr_Press(:,:,0)

    xy_SurfDSatQvapDTemp = EL * xy_SurfSatQvap / ( RVap * xy_SurfTemp * xy_SurfTemp )
    
    !-----------------------------------------------------------------
    !  フラックスの計算
    !  Calculate fluxes
    !-----------------------------------------------------------------
    xyr_UFlux(:,:,0) =   xyr_UFlux(:,:,0) - xy_SurfVelTransCoeff * xyz_U(:,:,0)

    xyr_VFlux(:,:,0) =   xyr_VFlux(:,:,0) - xy_SurfVelTransCoeff * xyz_V(:,:,0)

    xyr_TempFlux(:,:,0) =   xyr_TempFlux(:,:,0) + Cp * xy_SurfTempTransCoeff * (   xy_SurfTemp - xyz_Temp(:,:,0) / xy_SurfExner )

    xyr_QvapFlux(:,:,0) =   xyr_QvapFlux(:,:,0) + EL * xy_SurfQvapTransCoeff * xy_SurfHumidCoeff * ( xy_SurfSatQvap - xyz_Qvap(:,:,0) )
    
    !-----------------------------------------------------------------
    !  陰解行列の計算
    !  Calculate implicit matrices
    !-----------------------------------------------------------------
    xyaa_SurfTempMatrix = 0.0_DP
    xyaa_SurfQvapMatrix = 0.0_DP
    
    xy_SurfUVMatrix = xy_SurfVelTransCoeff
    
    xyaa_SurfTempMatrix(:,:,1,0) =   Cp * xy_SurfTempTransCoeff / xy_SurfExner
    xyaa_SurfTempMatrix(:,:,0,1) = - Cp * xy_SurfTempTransCoeff / xy_SurfExner
    
    xyaa_SurfQvapMatrix(:,:,1,0) =   Cp * xy_SurfQvapTransCoeff * xy_SurfHumidCoeff
    xyaa_SurfQvapMatrix(:,:,0,1) = - Cp * xy_SurfQvapTransCoeff * xy_SurfHumidCoeff
    
    do i = 0, imax-1
      do j = 0, jmax-1
        if ( xy_SurfCondition(i,j) >= 1 ) then
          
          xyaa_SurfTempMatrix(i,j,1,-1) = - Cp * xy_SurfTempTransCoeff(i,j)
          xyaa_SurfTempMatrix(i,j,0,0)  =   Cp * xy_SurfTempTransCoeff(i,j)
          
          xyaa_SurfQvapMatrix(i,j,1,-1) = - EL * xy_SurfQvapTransCoeff(i,j) * xy_SurfHumidCoeff(i,j) * xy_SurfDSatQvapDTemp(i,j)

          xyaa_SurfQvapMatrix(i,j,0,0)  = EL * xy_SurfQvapTransCoeff(i,j) * xy_SurfHumidCoeff(i,j) * xy_SurfDSatQvapDTemp(i,j)
          
        end if
      end do
    end do
    
    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhySurfFluxSurfaceFlux