subroutine PhyVerDiffClose( phy_vdif, err )
    !
    ! PHYVDIF 型の変数の終了処理を行います.
    ! なお, 与えられた *phy_vdif* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Deconstructor of "PHYVDIF".
    ! Note that if *phy_vdif* 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(PHYVDIF), intent(inout):: phy_vdif
    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 = 'PhyVerDiffClose'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

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

    !-----------------------------------------------------------------
    !  "PHYVDIF" の設定の消去
    !  Clear the settings for "PHYVDIF"
    !-----------------------------------------------------------------


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

  subroutine PhyVerDiffCreate( phy_vdif, imax, jmax, kmax, RAir, Cp, Grav, EL, FKarm, nmlfile, err )
    !
    ! PHYVDIF 型の変数の初期設定を行います.
    ! 他のサブルーチンを使用する前に必ずこのサブルーチンによって
    ! PHYVDIF 型の変数を初期設定してください.
    !
    ! なお, 与えられた *phy_vdif* が既に初期設定されている場合,
    ! プログラムはエラーを発生させます.
    !
    ! NAMELIST を利用する場合には引数 *nmlfile* に NAMELIST ファイル名
    ! を与えてください. NAMELIST 変数群の詳細に関しては 
    ! NAMELIST#phy_verdiff_nml を参照してください. 
    !
    ! Constructor of "PHYVDIF".
    ! Initialize *phy_vdif* by this subroutine, 
    ! before other procedures are used, 
    !
    ! Note that if *phy_vdif* is already initialized 
    ! by this procedure, error is occurred.
    !
    ! In order to use NAMELIST, specify a NAMELIST filename to 
    ! argument *nmlfile*. See "NAMELIST#phy_verdiff_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(PHYVDIF), intent(inout):: phy_vdif
    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):: RAir       ! $ R $ .      大気気体定数.   Gas constant of air
    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):: EL        ! $ L $ .      水の凝結の潜熱. Latent heat of condensation of water vapor
    real(DP), intent(in):: FKarm     ! $ k $ .                 カルマン定数.   Karman constant
    character(*), intent(in), optional:: nmlfile
                              ! NAMELIST ファイルの名称. 
                              ! この引数に空文字以外を与えた場合, 
                              ! 指定されたファイルから 
                              ! NAMELIST 変数群を読み込みます. 
                              ! ファイルを読み込めない場合にはエラーを
                              ! 生じます.
                              !
                              ! NAMELIST 変数群の詳細に関しては 
                              ! NAMELIST#phy_verdiff_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_verdiff_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 = 'PhyVerDiffCreate'
  continue
    call BeginSub( subname, version )
    stat = DC_NOERR
    cause_c = ''

    !-----------------------------------------------------------------
    !  初期設定のチェック
    !  Check initialization
    !-----------------------------------------------------------------
    if ( phy_vdif % initialized ) then
      stat = DC_EALREADYINIT
      cause_c = 'PHYVDIF'
      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

    !-----------------------------------------------------------------
    !  "PHYVDIF" の設定
    !  Configure the settings for "PHYVDIF"
    !-----------------------------------------------------------------
    phy_vdif % imax  = imax  
    phy_vdif % jmax  = jmax  
    phy_vdif % kmax  = kmax  
    phy_vdif % RAir  = RAir  
    phy_vdif % Cp    = Cp    
    phy_vdif % Grav  = Grav  
    phy_vdif % EL    = EL    
    phy_vdif % FKarm = FKarm 

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

    !-------------------------
    !  オプショナル引数からの値
    !  Values from optional arguments
!!$    phy_vdif % param_i = param_i
!!$    if ( present(param_r) )  phy_vdif % param_r = param_r
!!$    if ( present(param_c) )  phy_vdif % 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_vdif % param_i, &   ! (inout)
!!$        & param_r = phy_vdif % param_r, &   ! (inout)
!!$        & param_c_ = phy_vdif % 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_vdif % param_i < 0 ) then
!!$      stat = DC_ENEGATIVE
!!$      cause_c = 'param_i'
!!$      goto 999
!!$    end if


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

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

  subroutine PhyVerDiffNmlRead( 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_verdiff_nml/ &
!!$      & param_i, param_r, param_c
                              ! phy_verdiff モジュール用
                              ! NAMELIST 変数群名.
                              !
                              ! phy_verdiff#Create を使用する際に, 
                              ! オプショナル引数 *nmlfile* へ NAMELIST 
                              ! ファイル名を指定することで, そのファイルから
                              ! この NAMELIST 変数群を読み込みます.
                              !
                              ! NAMELIST group name for 
                              ! "phy_verdiff" module.
                              ! 
                              ! If a NAMELIST filename is specified to 
                              ! an optional argument *nmlfile* 
                              ! when "phy_verdiff#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 = 'PhyVerDiffNmlRead'
  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_verdiff_nml, iostat = iostat_nml ) ! (out)
!!$    if ( iostat_nml == 0 ) then
!!$      call MessageNotify( 'M', subname, &
!!$        & 'NAMELIST group "%c" is loaded from "%c".', &
!!$        & c1 = 'phy_verdiff_nml', c2 = trim(nmlfile) )
!!$      write(STDOUT, nml = phy_verdiff_nml)
!!$    else
!!$      call MessageNotify( 'W', subname, &
!!$        & 'NAMELIST group "%c" is not found in "%c" (iostat=%d).', &
!!$        & c1 = 'phy_verdiff_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 PhyVerDiffNmlRead

  subroutine PhyVerDiffPutLine( phy_vdif, unit, indent, err )
    !
    ! 引数 *phy_vdif* に設定されている情報を印字します.
    ! デフォルトではメッセージは標準出力に出力されます. 
    ! *unit* に装置番号を指定することで, 出力先を変更することが可能です.
    !
    ! Print information of *phy_vdif*.
    ! 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(PHYVDIF), intent(in):: phy_vdif
    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 = 'PhyVerDiffPutLine'
  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


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

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

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

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

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

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

  subroutine PhyVerDiffVerDiffCoefficients( phy_vdif, xyr_GeoPot, xyr_DVelDz, xyr_BulkRiNum, xyr_VelDiffCoeff, xyr_TempDiffCoeff, xyr_QvapDiffCoeff, err )
    !
    ! 拡散係数を計算します. 
    !
    ! なお, 与えられた *phy_vdif* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! Calculate diffusion coefficients. 
    !
    ! If *phy_vdif* 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(PHYVDIF), intent(inout):: phy_vdif
    real(DP), intent(in):: xyr_GeoPot (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! $ \phi $ . ジオポテンシャル (半整数レベル). 
                              ! Geo-potential (half level)
    real(DP), intent(in):: xyr_DVelDz (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! $ \DD{|v|}{z} $
    real(DP), intent(in):: xyr_BulkRiNum (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $
    real(DP), intent(out):: xyr_VelDiffCoeff (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP), intent(out):: xyr_TempDiffCoeff (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP), intent(out):: xyr_QvapDiffCoeff (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion 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. 


    !-----------------------------------
    !  格子点, 物理定数
    !  Grid number, physical constants
    integer:: imax ! 経度格子点数. 
                   ! Number of grid points in longitude
    integer:: jmax ! 緯度格子点数. 
                   ! Number of grid points in latitude
    integer:: kmax ! 鉛直層数. 
                   ! Number of vertical level
    real(DP):: FKarm     ! $ k $ .                 カルマン定数.   Karman constant

    !-----------------------------------
    !  拡散係数算出のための定数
    !  Constants for calculation of diffusion coefficients
    real(DP), parameter:: MixLengthMax      = 300.0_DP
                              ! 最大混合距離. 
                              ! Maximum mixing length
    real(DP), parameter:: TildeShMin        =   0.0_DP
                              ! $ \tilde{S_h} $ 最小値. 
                              ! Minimum $ \tilde{S_h} $
    real(DP), parameter:: TildeSmMin        =   0.0_DP
                              ! $ \tilde{S_m} $ 最小値. 
                              ! Minimum $ \tilde{S_m} $
    real(DP), parameter:: VelDiffCoeffMin  =   0.1_DP
                              ! $ \Dvect{u} $ 拡散係数最小値. 
                              ! Minimum diffusion coefficient of $ \Dvect{u} $
    real(DP), parameter:: TempDiffCoeffMin =   0.1_DP
                              ! $ T $ 拡散係数最小値. 
                              ! Minimum diffusion coefficient of $ T $
    real(DP), parameter:: QvapDiffCoeffMin =   0.1_DP
                              ! $ q $ 拡散係数最小値. 
                              ! Minimum diffusion coefficient of $ q $
    real(DP), parameter:: VelDiffCoeffMax  = 10000.0_DP
                              ! $ \Dvect{u} $ 拡散係数最大値. 
                              ! Maximum diffusion coefficient of $ \Dvect{u} $
    real(DP), parameter:: TempDiffCoeffMax = 10000.0_DP
                              ! $ T $ 拡散係数最大値. 
                              ! Maximum diffusion coefficient of $ T $
    real(DP), parameter:: QvapDiffCoeffMax = 10000.0_DP
                              ! $ q $ 拡散係数最大値. 
                              ! Maximum diffusion coefficient of $ q $

    real(DP):: xyr_FluxRiNum (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! フラックス $ R_i $ 数. 
                              ! Flux $ R_i $ number
    real(DP):: xyr_TildeSh (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! $ \tilde{S_h} $ (温度, 比湿). 
                              ! $ \tilde{S_h} $ (temperature, specific humidity)
    real(DP):: xyr_TildeSm (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! $ \tilde{S_m} $ (運動量). 
                              ! $ \tilde{S_m} $ (momentum)
    real(DP):: xyr_MixLength (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! 混合距離. 
                              ! Mixing length

    !-----------------------------------
    !  Mellor Yamada Level 2 定数
    !  Constants for Mellor Yamada Level 2
    real(DP), parameter:: ParamA1 = 0.92_DP
    real(DP), parameter:: ParamB1 = 16.6_DP
    real(DP), parameter:: ParamA2 = 0.74_DP
    real(DP), parameter:: ParamB2 = 10.1_DP
    real(DP), parameter:: ParamC1 = 0.08_DP

    real(DP):: Alpha1, Alpha2
    real(DP):: Beta1, Beta2, Beta3, Beta4
    real(DP):: Gamma1, Gamma2
    real(DP):: CriticalFluxRiNum

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: i, j, k         ! DO ループ用作業変数
                              ! Work variables for DO loop

    integer:: stat
    character(STRING):: cause_c
    character(*), parameter:: subname = 'PhyVerDiffVerDiffCoefficients'
  continue
    call BeginSub( subname )
    stat = DC_NOERR
    cause_c = ''

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

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


    !-----------------------------------------------------------------
    !  定数計算
    !  Calculate constants
    !-----------------------------------------------------------------
    Gamma1 = ( 1.0_DP / 3.0_DP ) - ( 2.0_DP * ParamA1 / ParamB1 )
    Gamma2 = ( ParamB2 / ParamB1 ) + ( 6.0_DP * ParamA1 / ParamB1 )
    Alpha1 = 3.0_DP  * ParamA2 * Gamma1
    Alpha2 = 3.0_DP  * ParamA2 * ( Gamma1 + Gamma1 )
    Beta1  = ParamA1 * ParamB1 * ( Gamma1 - ParamC1 )
    Beta2  = ParamA1 * (   ParamB1 * ( Gamma1 - ParamC1 ) + 6.0_DP * ParamA1 + 3.0_DP * ParamA2 )
    Beta3  = ParamA2 * ParamB1 * Gamma1
    Beta4  = ParamA2 * ( ParamB1 * ( Gamma1 + Gamma2 ) - 3.0_DP * ParamA1 )
    CriticalFluxRiNum = Gamma1 / ( Gamma1 + Gamma2 )
    
    !-----------------------------------------------------------------
    !  フラックス $ R_i $ 数の算出
    !  Calculate flux $ R_i $ number
    !-----------------------------------------------------------------
    xyr_FluxRiNum = (   Beta1 + Beta4 * xyr_BulkRiNum - sqrt(   ( Beta1 + Beta4 * xyr_BulkRiNum )**2 - 4.0_DP * Beta2 * Beta3 * xyr_BulkRiNum ) ) / ( 2.0_DP * Beta2 )
    
    !-----------------------------------------------------------------
    !  $ \tilde{S_h} $ と $ \tilde{S_m} $ の算出
    !  Calculate $ \tilde{S_h} $ and $ \tilde{S_m} $
    !-----------------------------------------------------------------
    xyr_TildeSh = 0.0_DP
    xyr_TildeSm = 0.0_DP
    
    do k = 0, kmax-1
      do i = 0, imax-1
        do j = 0, jmax-1
          
          if ( xyr_FluxRiNum(i,j,k) < CriticalFluxRiNum ) then 
            
            xyr_TildeSh(i,j,k) = (   Alpha1 - Alpha2 * xyr_FluxRiNum(i,j,k) ) / ( 1.0_DP - 1.0_DP * xyr_FluxRiNum(i,j,k) )

            xyr_TildeSm(i,j,k) = (   Beta1 - Beta2 * xyr_FluxRiNum(i,j,k) ) / ( Beta3 - Beta4 * xyr_FluxRiNum(i,j,k) ) * xyr_TildeSh(i,j,k)

            xyr_TildeSh(i,j,k) = max( xyr_TildeSh(i,j,k), TildeShMin )
            xyr_TildeSm(i,j,k) = max( xyr_TildeSm(i,j,k), TildeSmMin )
            
          else
            
            xyr_TildeSh(i,j,k) = TildeShMin
            xyr_TildeSm(i,j,k) = TildeSmMin
            
          end if
          
        end do
      end do
    end do
    
    
    !-----------------------------------------------------------------
    !  混合距離の算出
    !  Calculate mixing length
    !-----------------------------------------------------------------
    xyr_MixLength = FKarm * xyr_GeoPot / (1.0_DP + FKarm * xyr_GeoPot / MixLengthMax )
    
    !-----------------------------------------------------------------
    !  拡散係数の算出
    !  Calculate diffusion constants
    !-----------------------------------------------------------------
    
    xyr_VelDiffCoeff = xyr_MixLength**2 * xyr_DVelDz * sqrt ( ParamB1 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_TildeSm ) * xyr_TildeSm
    
    
    xyr_TempDiffCoeff = xyr_MixLength ** 2 * xyr_DVelDz * sqrt ( ParamB1 * ( 1.0_DP - xyr_FluxRiNum ) * xyr_TildeSm ) * xyr_TildeSh
    
    xyr_QvapDiffCoeff = xyr_TempDiffCoeff
    
    do k = 0, kmax-1
      do i = 0, imax-1
        do j = 0, jmax-1
          xyr_VelDiffCoeff(i,j,k) = max( min( xyr_VelDiffCoeff(i,j,k), VelDiffCoeffMax ), VelDiffCoeffMin )
          xyr_TempDiffCoeff(i,j,k) = max( min( xyr_TempDiffCoeff(i,j,k), TempDiffCoeffMax ), TempDiffCoeffMin )
          xyr_QvapDiffCoeff(i,j,k) = max( min( xyr_QvapDiffCoeff(i,j,k), QvapDiffCoeffMax ), QvapDiffCoeffMin )
        end do
      end do
    end do
    
    xyr_VelDiffCoeff(:,:,0)     = 0.0_DP
    xyr_TempDiffCoeff(:,:,0)    = 0.0_DP
    xyr_QvapDiffCoeff(:,:,0)    = 0.0_DP
    xyr_VelDiffCoeff(:,:,kmax)  = 0.0_DP
    xyr_TempDiffCoeff(:,:,kmax) = 0.0_DP
    xyr_QvapDiffCoeff(:,:,kmax) = 0.0_DP

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

  subroutine PhyVerDiffVerticalDiffusion( phy_vdif, xyz_U, xyz_V, xyz_Temp, xyr_Temp, xyz_QVap, xyz_Press, xyr_Press, xyz_GeoPot, xyr_GeoPot, xyr_UFlux, xyr_VFlux, xyr_TempFlux, xyr_QVapFlux, xyza_UVMatrix, xyza_TempMatrix, xyza_QVapMatrix, err )
    !
    ! なお, 与えられた *phy_vdif* が Create によって初期設定
    ! されていない場合, プログラムはエラーを発生させます.
    !
    ! If *phy_vdif* 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(PHYVDIF), intent(inout):: phy_vdif
    real(DP), intent(in):: xyz_U (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax-1)
                              ! $ U $ . 東西風速. Zonal wind
    real(DP), intent(in):: xyz_V (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax-1)
                              ! $ V $ . 南北風速. Meridional wind

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

!!$    integer:: param_i
!!$    real(DP):: param_r
!!$    character(STRING):: param_c

    !-----------------------------------
    !  作業変数
    !  Work variables
    integer:: kmax ! 鉛直層数. 
                   ! Number of vertical level
    real(DP):: RAir       ! $ R $ .      大気気体定数.   Gas constant of air
    real(DP):: Grav       ! $ g $ .      重力加速度.     Gravitational acceleration
    real(DP):: Cp        ! $ C_p $ .    大気定圧比熱.   Specific heat of air at constant pressure
    real(DP):: EL        ! $ L $ .      水の凝結の潜熱. Latent heat of condensation of water vapor

    real(DP), parameter:: RefPress        = 1.0e+6_DP
                              ! 参照気圧. 
                              ! Reference air pressure
    real(DP), parameter:: BasePotTemp     = 300.0_DP
                              ! 基本温位. 
                              ! Base potential temperature
    real(DP), parameter:: SquareVelMin    = 0.1_DP
                              ! 風二乗差最小値. 
                              ! Minimum value of square of velocity
    real(DP), parameter:: BulkRiNumMin    = - 100.0_DP
                              ! バルク $ R_i $ 数最小値. 
                              ! Minimum value of bulk $ R_i $

    real(DP):: xyr_DVelDz (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! $ \DD{|\Dvect{v}|}{z} $
    real(DP):: xyr_BulkRiNum (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! バルク $ R_i $ 数. 
                              ! Bulk $ R_i $
    real(DP):: xyr_VelTransCoeff (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! 輸送係数：運動量. 
                              ! Transfer coefficient: velocity
    real(DP):: xyr_TempTransCoeff (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! 輸送係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP):: xyr_QvapTransCoeff (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! 輸送係数：比湿. 
                              ! Transfer coefficient: specific humidity
    real(DP):: xyr_VelDiffCoeff (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! 拡散係数：運動量. 
                              ! Diffusion coefficient: velocity
    real(DP):: xyr_TempDiffCoeff (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! 拡散係数：温度. 
                              ! Transfer coefficient: temperature
    real(DP):: xyr_QvapDiffCoeff (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! 拡散係数：比湿. 
                              ! Diffusion coefficient: specific humidity
    real(DP):: xyz_Exner (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax-1)
                              ! Exner 関数 (整数レベル). 
                              ! Exner function (full level)
    real(DP):: xyr_Exner (0:phy_vdif%imax-1, 0:phy_vdif%jmax-1, 0:phy_vdif%kmax)
                              ! Exner 関数 (半整数レベル). 
                              ! Exner function (half level)


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

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

    !-----------------------------------------------------------------
    !  *phy_vdif* に格納されている設定値の取り出し
    !  Fetch setting values stored in *phy_vdif*
    !-----------------------------------------------------------------
    kmax  = phy_vdif % kmax  
    RAir  = phy_vdif % RAir  
    Cp    = phy_vdif % Cp    
    Grav  = phy_vdif % Grav  
    EL    = phy_vdif % EL    



    !-----------------------------------------------------------------
    !  Exner 関数算出
    !  Calculate Exner functions
    !-----------------------------------------------------------------
    xyz_Exner = ( xyz_Press / RefPress ) ** ( RAir / Cp )
    xyr_Exner = ( xyr_Press / RefPress ) ** ( RAir / Cp )

    !-----------------------------------------------------------------
    !  バルク $ R_i $ 数算出
    !  Calculate bulk $ R_i $
    !-----------------------------------------------------------------
    xyr_DVelDz    = 0.0_DP
    xyr_BulkRiNum = 0.0_DP
    
    do k = 1, kmax-1
      xyr_DVelDz(:,:,k) = sqrt( max( SquareVelMin , ( xyz_U(:,:,k) - xyz_U(:,:,k-1) )**2 + ( xyz_V(:,:,k) - xyz_V(:,:,k-1) )**2 ) ) / ( xyz_GeoPot(:,:,k) - xyz_GeoPot(:,:,k-1) )
      
      xyr_BulkRiNum(:,:,k) = Grav / BasePotTemp * (   xyz_Temp(:,:,k)   / xyz_Exner(:,:,k) - xyz_Temp(:,:,k-1) / xyz_Exner(:,:,k-1) ) / ( xyz_GeoPot(:,:,k) - xyz_GeoPot(:,:,k-1) ) / xyr_DVelDz(:,:,k)**2
      
      xyr_BulkRiNum(:,:,k) = max( xyr_BulkRiNum(:,:,k) , BulkRiNumMin )
      
    end do
    
    !-----------------------------------------------------------------
    !  拡散係数の計算
    !  Calculate diffusion coefficients
    !-----------------------------------------------------------------
    call VerDiffCoefficients( phy_vdif = phy_vdif, xyr_GeoPot = xyr_GeoPot, xyr_DVelDz = xyr_DVelDz, xyr_BulkRiNum = xyr_BulkRiNum, xyr_VelDiffCoeff = xyr_VelDiffCoeff, xyr_TempDiffCoeff = xyr_TempDiffCoeff, xyr_QvapDiffCoeff = xyr_QvapDiffCoeff )      ! (out)
    
    !-----------------------------------------------------------------
    !  浅い積雲対流
    !  Shallow cumulus convection
    !-----------------------------------------------------------------
    
    
    !-----------------------------------------------------------------
    !  拡散係数の出力
    !  Output diffusion coefficients
    !-----------------------------------------------------------------
    
    
    !-----------------------------------------------------------------
    !  輸送係数の計算
    !  Calculate transfer coefficient
    !-----------------------------------------------------------------
    xyr_VelTransCoeff  = 0.0_DP
    xyr_TempTransCoeff = 0.0_DP
    xyr_QvapTransCoeff = 0.0_DP
    
    do k = 1, kmax-1
      xyr_VelTransCoeff(:,:,k) = xyr_VelDiffCoeff(:,:,k) * xyr_Press(:,:,k) / RAir / xyr_Temp(:,:,k) / ( xyz_GeoPot(:,:,k) - xyz_GeoPot(:,:,k-1) )

      xyr_TempTransCoeff(:,:,k) = xyr_TempDiffCoeff(:,:,k) * xyr_Press(:,:,k) / RAir / xyr_Temp(:,:,k) / ( xyz_GeoPot(:,:,k) - xyz_GeoPot(:,:,k-1) )

      xyr_QvapTransCoeff(:,:,k) = xyr_QvapDiffCoeff(:,:,k) * xyr_Press(:,:,k) / RAir / xyr_Temp(:,:,k) / ( xyz_GeoPot(:,:,k) - xyz_GeoPot(:,:,k-1) )
    end do
    
    !-----------------------------------------------------------------
    !  フラックスの計算
    !  Calculate fluxes
    !-----------------------------------------------------------------
    do k = 1, kmax-1
      xyr_UFlux(:,:,k) = xyr_UFlux(:,:,k) + xyr_VelTransCoeff(:,:,k) * ( xyz_U(:,:,k-1) - xyz_U(:,:,k) )
      
      xyr_VFlux(:,:,k) = xyr_VFlux(:,:,k) + xyr_VelTransCoeff(:,:,k) * ( xyz_V(:,:,k-1) - xyz_V(:,:,k) )
      
      xyr_TempFlux(:,:,k) = xyr_TempFlux(:,:,k) + Cp * xyr_TempTransCoeff(:,:,k) * xyr_Exner(:,:,k) * (   xyz_Temp(:,:,k-1) / xyz_Exner(:,:,k-1) - xyz_Temp(:,:,k)   / xyz_Exner(:,:,k)  )
      
      xyr_QvapFlux(:,:,k) = xyr_QvapFlux(:,:,k) + EL * xyr_QvapTransCoeff(:,:,k) * ( xyz_Qvap(:,:,k-1) - xyz_Qvap(:,:,k) )
      
    end do
    
    !-----------------------------------------------------------------
    !  陰解行列の計算
    !  Calculate implicit matrices
    !-----------------------------------------------------------------
    do k = 1, kmax-1
      xyza_UVMatrix(:,:,k,0)  = xyza_UVMatrix(:,:,k,0) + xyr_VelTransCoeff(:,:,k)

      xyza_UVMatrix(:,:,k,-1) = - xyr_VelTransCoeff(:,:,k)
      
      xyza_TempMatrix(:,:,k,0)  = xyza_TempMatrix(:,:,k,0) + Cp * xyr_TempTransCoeff(:,:,k) * xyr_Exner(:,:,k) / xyz_Exner(:,:,k)

      xyza_TempMatrix(:,:,k,-1) = - Cp * xyr_TempTransCoeff(:,:,k) * xyr_Exner(:,:,k) / xyz_Exner(:,:,k-1)
      
      xyza_QvapMatrix(:,:,k,0) = xyza_QvapMatrix(:,:,k,0) + Cp * xyr_QvapTransCoeff(:,:,k)

      xyza_QvapMatrix(:,:,k,-1) = - Cp * xyr_QvapTransCoeff(:,:,k)
    end do
    
    do k = 0, kmax-2
      xyza_UVMatrix(:,:,k,0) = xyza_UVMatrix(:,:,k,0) + xyr_VelTransCoeff(:,:,k+1)

      xyza_UVMatrix(:,:,k,1) = - xyr_VelTransCoeff(:,:,k+1)
      
      xyza_TempMatrix(:,:,k,0)  = xyza_TempMatrix(:,:,k,0) + Cp * xyr_TempTransCoeff(:,:,k+1) * xyr_Exner(:,:,k+1) / xyz_Exner(:,:,k)

      xyza_TempMatrix(:,:,k,1)  = - Cp * xyr_TempTransCoeff(:,:,k+1) * xyr_Exner(:,:,k+1) / xyz_Exner(:,:,k+1)
      
      xyza_QvapMatrix(:,:,k,0) = xyza_QvapMatrix(:,:,k,0) + Cp * xyr_QvapTransCoeff(:,:,k+1)

      xyza_QvapMatrix(:,:,k,1) = - Cp * xyr_QvapTransCoeff(:,:,k+1)
    end do
    
    !-----------------------------------------------------------------
    !  終了処理, 例外処理
    !  Termination and Exception handling
    !-----------------------------------------------------------------
999 continue
    call StoreError( stat, subname, err, cause_c )
    call EndSub( subname )
  end subroutine PhyVerDiffVerticalDiffusion