subroutine VenusSimpleForcing( xyz_UB, xyz_VB, xyz_TempB, xy_PsB, xyz_Press, xyr_Press, xyr_Temp, xyz_Height, xyz_Exner, xyr_Exner, xyz_DUDt, xyz_DVDt, xyz_DTempDt )

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

    real(DP), intent(in ) :: xyz_UB     (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_VB     (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_TempB  (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xy_PsB     (0:imax-1,1:jmax)
    real(DP), intent(in ) :: xyz_Press  (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyr_Press  (0:imax-1,1:jmax,0:kmax)
    real(DP), intent(in ) :: xyr_Temp   (0:imax-1,1:jmax,0:kmax)
    real(DP), intent(in ) :: xyz_Height (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_Exner  (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyr_Exner  (0:imax-1,1:jmax,0:kmax)
    real(DP), intent(out) :: xyz_DUDt   (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_DVDt   (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_DTempDt(0:imax-1,1:jmax,1:kmax)


    !
    ! local variables
    !
    real(DP) :: xyz_DUDtSFCFric (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DVDtSFCFric (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DUDtVDiff   (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DVDtVDiff   (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DTempDtVDiff(0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DTempDtRadL (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DTempDtRadS (0:imax-1,1:jmax,1:kmax)


    ! 初期化
    ! Initialization
    !
    if ( .not. venus_simple_forcing_inited ) call VenusSimpleForcingInit

    call VenusSimpleRadForcing( xy_PsB, xyz_Press, xyz_TempB, xyz_Height, xyz_DTempDtRadL, xyz_DTempDtRadS )


    call VenusSimpleSurfFriction( xyz_UB, xyz_VB, xyz_DUDtSFCFric, xyz_DVDtSFCFric )

    call VenusSimpleVDiff( xyz_UB, xyz_VB, xyz_TempB, xyr_Press, xyr_Temp, xyz_Height, xyz_Exner, xyr_Exner, xyz_DUDtVDiff, xyz_DVDtVDiff, xyz_DTempDtVDiff )


    xyz_DUDt    = xyz_DUDtVDiff    + xyz_DUDtSFCFric
    xyz_DVDt    = xyz_DVDtVDiff    + xyz_DVDtSFCFric

    xyz_DTempDt = xyz_DTempDtVDiff + xyz_DTempDtRadL + xyz_DTempDtRadS


  end subroutine VenusSimpleForcing

  subroutine VenusSimpleDTempDtRadS( y_CosLat, xyz_Press, xyz_Height, xyz_DTempDtRadS )

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, GasRDry

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: DP, STRING     ! 文字列.       Strings.


    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

    real(DP), intent(in ) :: y_CosLat(1:jmax)
    real(DP), intent(in ) :: xyz_Press      (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_Height     (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_DTempDtRadS(0:imax-1,1:jmax,1:kmax)


    !
    ! local variables
    !
    real(DP)   :: scaleheight
    real(DP)   :: DTempDtRadSMax
    integer(4) :: i, j, k


!!$    xyz_DTempDtRadS &
!!$      & = 5.0d0 / dayearth * exp( - ( ( xyz_Height - 55.0d3 ) / 10.0d3 )**2  )
!!$
!!$    do k = 1, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          if( xyz_Height(i,j,k) .le. 55.0d3 ) then
!!$            if( xyz_DTempDtRadS(i,j,k) .lt. 0.5d0 / dayearth ) then
!!$              xyz_DTempDtRadS(i,j,k) = 0.5d0 / dayearth
!!$            end if
!!$          end if
!!$        end do
!!$      end do
!!$    end do


    scaleheight = GasRDry * 300.0d0 / Grav

    xyz_DTempDtRadS = 5.0d0 / dayearth * exp( - ( ( - scaleheight * log( xyz_Press / 500.0d2 ) ) / ( 2.0d0 * scaleheight ) )**2 )

    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_Press(i,j,k) > 500.0d2 ) then
            if ( xyz_DTempDtRadS(i,j,k) .lt. 0.5d0 / dayearth ) then
              xyz_DTempDtRadS(i,j,k) = 0.5d0 / dayearth
            end if
          end if
        end do
      end do
    end do




!!$    do k = 1, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1


!!$          if( xyz_Press(i,j,k) .le. 1.0d5 ) then
!!$!                  gswrh( i, j, k ) = 5.0d0 / dayearth
!!$            xyz_DTempDtRadS(i,j,k) = 5.0d0 / dayearth &
!!$              & * exp( - ( 5.0d3 * log( xyz_Press(i,j,k) / 1.0d5 ) / 15.0d3 )**2  )
!!$          else
!!$            xyz_DTempDtRadS(i,j,k) &
!!$              & = log( ( 5.0d0 / dayearth ) / ( 1.0d-4 / dayearth ) ) &
!!$              & / log(   1.0d5              /   100.0d5             ) &
!!$              & * log(   xyz_Press(i,j,k)   /   100.0d5             ) &
!!$              & + log(   1.0d-4 / dayearth  )
!!$            xyz_DTempDtRadS(i,j,k) = exp( xyz_DTempDtRadS(i,j,k) )
!!$            if( xyz_DTempDtRadS(i,j,k) .lt. 0.5d0 / dayearth ) then
!!$              xyz_DTempDtRadS(i,j,k) = 0.5d0 / dayearth
!!$            end if
!!$          end if


          !-----


!!$          DTempDtRadSMax = 3.0d0 / dayearth
!!$
!!$          if( xyz_Press(i,j,k) .le. 1.0d4 ) then
!!$            xyz_DTempDtRadS(i,j,k) = DTempDtRadSMax &
!!$              & * exp( - ( 5.0d3 * log( xyz_Press(i,j,k) / 1.0d4 ) / 10.0d3 )**2  )
!!$          else if( xyz_Press(i,j,k) .le. 1.0d5 ) then
!!$            xyz_DTempDtRadS(i,j,k) = DTempDtRadSMax
!!$
!!$!               if( gp( i, j, k ) .le. 1.0d5 ) then
!!$!                  gswrh( i, j, k ) = sw_hr_peak &
!!$!                       * exp( - ( 5.0d3 * log( gp( i, j, k ) / 1.0d5 ) / 15.0d3 )**2  )
!!$
!!$          else
!!$            xyz_DTempDtRadS(i,j,k) &
!!$              & = log( DTempDtRadSMax       / ( 1.0d-4 / dayearth ) ) &
!!$              & / log(   1.0d5              /   100.0d5             ) &
!!$              & * log( xyz_Press(i,j,k)     /   100.0d5             ) &
!!$              & + log(   1.0d-4 / dayearth  )
!!$            xyz_DTempDtRadS(i,j,k) = exp( xyz_DTempDtRadS(i,j,k) )
!!$            if( xyz_DTempDtRadS(i,j,k) .lt. 0.5d0 / dayearth ) then
!!$              xyz_DTempDtRadS(i,j,k) = 0.5d0 / dayearth
!!$            end if
!!$!                  if( gswrh( i, j, k ) .lt. 0.15d0 / dayearth ) then
!!$!                     gswrh( i, j, k ) = 0.15d0 / dayearth
!!$!                  end if
!!$          end if


!!$        end do
!!$      end do
!!$    end do


    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          xyz_DTempDtRadS(i,j,k) = xyz_DTempDtRadS(i,j,k) * y_CosLat(j)
        end do
      end do
    end do


  end subroutine VenusSimpleDTempDtRadS

  subroutine VenusSimpleForcingInit
    !
    ! venus_simple_forcing モジュールの初期化を行います.
    ! NAMELIST#venus_simple_forcing_nml の読み込みはこの手続きで行われます.
    !
    ! "venus_simple_forcing" module is initialized.
    ! "NAMELIST#venus_simple_forcing_nml" is loaded in this procedure.
    !


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

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: GasRDry, CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: y_Lat, z_Sigma               ! $ \sigma $ レベル (整数).
                              ! Full $ \sigma $ level

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output

    ! 文字列操作
    ! Character handling
    !
    use dc_string, only: StoA

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable

    ! 宣言文 ; Declaration statements
    !
    implicit none

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.
                              ! IOSTAT of NAMELIST read

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /venus_simple_forcing_nml/ SurfFrictionTimeConstInEarthDay, VDiffCoefM, VDiffCoefH, FlagConstNCC, ConstNCCInEarthDay
          !
          ! デフォルト値については初期化手続 "venus_simple_forcing#VenusSimpleForcingInit"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "venus_simple_forcing#VenusSimpleForcingInit" for the default values.
          !

    ! 実行文 ; Executable statement
    !

    if ( venus_simple_forcing_inited ) return
!!$    call InitCheck


    ! デフォルト値の設定
    ! Default values settings
    !
    SurfFrictionTimeConstInEarthDay = 30.0d0
    VDiffCoefM                      =  1.0d0
    VDiffCoefH                      =  1.0d0
    FlagConstNCC                    = .true.
    ConstNCCInEarthDay              = 30.0d0


    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = venus_simple_forcing_nml, iostat = iostat_nml )
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
!$      if ( iostat_nml == 0 ) write( STDOUT, nml = venus_simple_forcing_nml )
    end if


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'VTempEq', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'radiative equilibrium temperature', 'K' )
    call HistoryAutoAddVariable( 'VSRadHR', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'solar heating rate', 'K s-1' )
    call HistoryAutoAddVariable( 'VEquivTempEq', (/ 'lon ', 'lat ', 'sig ', 'time' /), '"equivalent" radiative equilibrium temperature', 'K' )
    call HistoryAutoAddVariable( 'VUBalance', (/ 'lon ', 'lat ', 'sig ', 'time' /), 'balanced zonal wind', 'm s-1' )

    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  SurfFrictionTimeConstInEarthDay = %f', d = (/ SurfFrictionTimeConstInEarthDay /) )
    call MessageNotify( 'M', module_name, '  VDiffCoefM                      = %f', d = (/ VDiffCoefM /) )
    call MessageNotify( 'M', module_name, '  VDiffCoefH                      = %f', d = (/ VDiffCoefH /) )
    call MessageNotify( 'M', module_name, '  FlagConstNCC                    = %b', l = (/ FlagConstNCC /) )
    call MessageNotify( 'M', module_name, '  ConstNCCInEarthDay              = %f', d = (/ ConstNCCInEarthDay /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    venus_simple_forcing_inited = .true.


  end subroutine VenusSimpleForcingInit

  subroutine VenusSimpleNCCoef( xy_Ps, xyz_Press, xyz_NCC )

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: DP, STRING     ! 文字列.       Strings.


    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

    real(DP), intent(in ) :: xy_Ps    (0:imax-1,1:jmax)
    real(DP), intent(in ) :: xyz_Press(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_NCC  (0:imax-1,1:jmax,1:kmax)


    !
    ! local variables
    !
    real(DP) :: xyz_alp1  (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_alp2  (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_alp3  (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_lnPRat(0:imax-1,1:jmax,1:kmax)
    real(DP) :: NCTimeConst
    real(DP) :: NCTimeConst0
    integer(4) :: i, j, k


    if ( FlagConstNCC ) then

      xyz_NCC = 1.0d0 / ( ConstNCCInEarthDay * dayearth )

    else

      ! Thermal damping by Hou and Farrel (1987)
      !
      do k = 1, kmax
        xyz_lnPRat(:,:,k) = log( xyz_Press(:,:,k) / xy_Ps(:,:) )
      end do
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            if( -xyz_lnPRat(i,j,k) .le. 5.0d0 ) then
              xyz_alp1(i,j,k) =  0.0d0
              xyz_alp2(i,j,k) =  0.9d0
              xyz_alp3(i,j,k) =  0.0d0
            else if( -xyz_lnPRat(i,j,k) .le. 7.0d0 ) then
              xyz_alp1(i,j,k) = -4.5d0
              xyz_alp2(i,j,k) =  2.0d0
              xyz_alp3(i,j,k) =  5.0d0
            else
              xyz_alp1(i,j,k) = -8.5d0
              xyz_alp2(i,j,k) =  0.5d0
              xyz_alp3(i,j,k) =  7.0d0
            end if
          end do
        end do
      end do
      NCTimeConst0 = 1.32d9
      do k = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            NCTimeConst = NCTimeConst0 * exp( xyz_alp1(i,j,k) - xyz_alp2(i,j,k) * ( -xyz_lnPRat(i,j,k) - xyz_alp3(i,j,k) ) )
            xyz_NCC(i,j,k) = 1.0d0 / NCTimeConst

            if( xyz_NCC(i,j,k) .lt. 1.0d0 / ( 30.0d0 * dayearth ) ) then
              xyz_NCC(i,j,k) = 1.0d0 / ( 30.0d0 * dayearth )
            end if
          end do
        end do
      end do

    end if


  end subroutine VenusSimpleNCCoef

  subroutine VenusSimpleNCTempEq( xyz_Height, xyz_TempEq )

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: DP, STRING     ! 文字列.       Strings.


    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure

    use axesset, only: y_Lat, z_Sigma               ! $ \sigma $ レベル (整数).
                              ! Full $ \sigma $ level

    real(DP), intent(in ) :: xyz_Height(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_TempEq(0:imax-1,1:jmax,1:kmax )


    !
    ! local variables
    !
    real(DP)   :: SurfTemp
    real(DP)   :: z( 5 ), a( 6 ), ah( 5 ), d( 5 )
    integer(4) :: l


    ! Coefficients for thermal structure by Hou and Farrel (1987)
    !
!!$    z ( 1 ) =   0.0d3
!!$    z ( 2 ) =  10.0d3
!!$    z ( 3 ) =  25.0d3
!!$    z ( 4 ) =  55.0d3
!!$    z ( 5 ) = 100.0d3
!!$
!!$    ah( 1 ) =  -1.0d-3
!!$    ah( 2 ) =  -1.0d-3
!!$    ah( 3 ) =  -3.1d-3
!!$    ah( 4 ) =  -6.75d-3
!!$    ah( 5 ) =  10.0d-3
!!$
!!$    d ( 1 ) =  10.0d3
!!$    d ( 2 ) =  10.0d3
!!$    d ( 3 ) =   8.0d3
!!$    d ( 4 ) =   5.0d3
!!$    d ( 5 ) =  70.0d3


    ! Slightly modified coefficients for thermal structure by Hou and Farrel (1987)
    !
    z ( 1 ) =   0.0d3
    z ( 2 ) =  10.0d3
    z ( 3 ) =  25.0d3
!!$    z ( 4 ) =  55.0d3
    z ( 4 ) =  50.0d3
    z ( 5 ) = 100.0d3

    ah( 1 ) =  -1.0d-3
    ah( 2 ) =  -1.0d-3
!!$    ah( 3 ) =  -3.1d-3
    ah( 3 ) =  -2.0d-3
!!$    ah( 4 ) =  -6.75d-3
    ah( 4 ) =  -3.0d-3
    ah( 5 ) =  10.0d-3

    d ( 1 ) =  10.0d3
    d ( 2 ) =  10.0d3
!!$    d ( 3 ) =   8.0d3
    d ( 3 ) =  15.0d3
!!$    d ( 4 ) =   5.0d3
    d ( 4 ) =  10.0d3
    d ( 5 ) =  70.0d3



    a ( 1 ) =   0.0d0

    do l = 2, 6
      a( l ) = 2.0d0 * ah( l-1 ) * d( l-1 ) + a( l-1 )
    end do


    SurfTemp = 750.0d0
    xyz_TempEq = SurfTemp - Grav / CpDry * xyz_Height

    do l = 1, 5
!!$      if ( l == 4 ) cycle
      xyz_TempEq = xyz_TempEq - ( a(l+1) - a(l) ) * 0.5d0 * ( 1.0d0 + tanh( ( 0.0d0      - z(l) ) / d(l) ) )
      xyz_TempEq = xyz_TempEq + ( a(l+1) - a(l) ) * 0.5d0 * ( 1.0d0 + tanh( ( xyz_Height - z(l) ) / d(l) ) )
    end do

!!$    do l = 1, kmax
!!$      write( 90, * ) xyz_TempEq(0,jmax/2+1,l), z_sigma(l)
!!$    end do
!!$    call flush( 90 )
!!$    stop


  end subroutine VenusSimpleNCTempEq

  subroutine VenusSimpleRadForcing( xy_Ps, xyz_Press, xyz_Temp, xyz_Height, xyz_DTempDtRadL, xyz_DTempDtRadS )

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: DelTime, TimeN, TimesetClockStart, TimesetClockStop

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, GasRDry
                              ! $ R $ [J kg-1 K-1].
                              ! 乾燥大気の気体定数.
                              ! Gas constant of air

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: y_Lat, z_Sigma               ! $ \sigma $ レベル (整数).
                              ! Full $ \sigma $ level


    real(DP), intent(in ):: xy_Ps          (0:imax-1,1:jmax)
    real(DP), intent(in ):: xyz_Press      (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ):: xyz_Temp       (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ):: xyz_Height     (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out):: xyz_DTempDtRadL(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out):: xyz_DTempDtRadS(0:imax-1,1:jmax,1:kmax)


    !
    ! local variables
    !
    real(DP) :: y_CosLat(1:jmax)
    real(DP) :: xyz_TempEq(0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_NCC   (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_EquivTempEq(0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_Geopot(0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_UBalance(0:imax-1,1:jmax,1:kmax)
    integer  :: i, j, k
    integer  :: jp, jn


    y_CosLat = cos( y_Lat )

    call VenusSimpleNCTempEq( xyz_Height, xyz_TempEq )

    call VenusSimpleNCCoef( xy_Ps, xyz_Press, xyz_NCC )

    xyz_DTempDtRadL = - xyz_NCC * ( xyz_Temp - xyz_TempEq )


    call VenusSimpleDTempDtRadS( y_CosLat, xyz_Press, xyz_Height, xyz_DTempDtRadS )



!!$    do k = 1, kmax
!!$      do j = 1, jmax
!!$        write( 60, * ) j, xyz_Press(1,j,k), xyz_DTempDtRadS(1,j,k) * dayearth
!!$      end do
!!$      write( 60, * )
!!$    end do
!!$    call flush( 60 )
!!$
!!$    do k = 1, kmax
!!$      write( 61, * ) k, xyz_Height(1,1,k) * 1.0d-3, xyz_Press(1,1,k), 1.0d0 / xyz_NCC(1,1,k) / dayearth, xyz_TempEq(1,1,k), xyz_DTempDtRadS(1,jmax/2,k) * dayearth
!!$    end do
!!$    stop

    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_NCC(i,j,k) /= 0.0d0 ) then
            xyz_EquivTempEq(i,j,k) = xyz_TempEq(i,j,k) + xyz_DTempDtRadS(i,j,k) / xyz_NCC(i,j,k)
          else
            xyz_EquivTempEq(i,j,k) = xyz_TempEq(i,j,k)
          end if
        end do
      end do
    end do


    ! dp/dz = -rho g
    ! dp / dphi = -rho
    ! dphi / dp = -1/rho = - R T / p
    ! p dphi / dp = -1/rho = - R T
    ! dphi / dlogp = - R T

    k = 1
    xyz_Geopot(:,:,k) = 0.0d0 - GasRDry * xyz_EquivTempEq(:,:,k) * log( z_Sigma(k) )
    do k = 2, kmax
      xyz_Geopot(:,:,k) = xyz_Geopot(:,:,k-1) - GasRDry * ( xyz_EquivTempEq(:,:,k-1) + xyz_EquivTempEq(:,:,k) ) * 0.5d0 * log( z_Sigma(k) / z_Sigma(k-1) )
    end do

!!$    xyz_UBalance = xyz_Geopot
!!$
    do k = 1, kmax
      do j = 1, jmax
        if ( j == 1 ) then
          jp = 1
          jn = j + 1
        else if ( j == jmax ) then
          jp = j - 1
          jn = jmax
        else
          jp = j - 1
          jn = j + 1
        end if
        xyz_UBalance(:,j,k) = sqrt( - ( xyz_Geopot(:,jn,k) - xyz_Geopot(:,jp,k) ) / ( y_Lat(jn)          - y_Lat(jp)          ) / tan( y_Lat(j) ) )
      end do
    end do


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'VTempEq'     , xyz_TempEq )
    call HistoryAutoPut( TimeN, 'VSRadHR'     , xyz_DTempDtRadS )
    call HistoryAutoPut( TimeN, 'VEquivTempEq', xyz_EquivTempEq )
    call HistoryAutoPut( TimeN, 'VUBalance'   , xyz_UBalance )


  end subroutine VenusSimpleRadForcing

  subroutine VenusSimpleRayleighFriction( xyz_Press, xyz_UB, xyz_VB, xyz_DUDt, xyz_DVDt )

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level


    real(DP), intent(in   ) :: xyz_Press(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in   ) :: xyz_UB   (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in   ) :: xyz_VB   (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(inout) :: xyz_DUDt (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(inout) :: xyz_DVDt (0:imax-1,1:jmax,1:kmax)


    !
    ! local variables
    !
    real(DP)  :: yz_ZMU(1:jmax,1:kmax)
    real(DP)  :: yz_ZMV(1:jmax,1:kmax)
    integer(4):: i, j, k


    yz_ZMU = 0.0d0
    yz_ZMV = 0.0d0
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          yz_ZMU(j,k) = yz_ZMU(j,k) + xyz_UB(i,j,k)
          yz_ZMV(j,k) = yz_ZMV(j,k) + xyz_VB(i,j,k)
        end do
      end do
    end do
    yz_ZMU = yz_ZMU / dble( imax )
    yz_ZMV = yz_ZMV / dble( imax )

    !
    ! Rayleigh friction
    !
    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          xyz_DUDt(i,j,k) = xyz_DUDt(i,j,k) - ( xyz_UB(i,j,k) - yz_ZMU(j,k) ) / ( dayearth * ( xyz_Press(i,j,k) / xyz_Press(i,j,kmax) )**2 )
          xyz_DVDt(i,j,k) = xyz_DVDt(i,j,k) - ( xyz_VB(i,j,k) - yz_ZMV(j,k) ) / ( dayearth * ( xyz_Press(i,j,k) / xyz_Press(i,j,kmax) )**2 )
        end do
      end do
    end do


  end subroutine VenusSimpleRayleighFriction

  subroutine VenusSimpleSurfFriction( xyz_UB, xyz_VB, xyz_DUDt, xyz_DVDt )

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure


    use axesset  , only : y_Lat


    real(DP), intent(in ):: xyz_UB  (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ):: xyz_VB  (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out):: xyz_DUDt(0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out):: xyz_DVDt(0:imax-1,1:jmax,1:kmax)


    !
    ! local variables
    !
    real(DP) :: SurfFrictionTimeConst


    SurfFrictionTimeConst = SurfFrictionTimeConstInEarthDay * dayearth

    xyz_DUDt(:,:,2:kmax) = 0.0d0
    xyz_DVDt(:,:,2:kmax) = 0.0d0

    xyz_DUDt(:,:,1) = - xyz_UB(:,:,1) / SurfFrictionTimeConst
    xyz_DVDt(:,:,1) = - xyz_VB(:,:,1) / SurfFrictionTimeConst


  end subroutine VenusSimpleSurfFriction

  subroutine VenusSimpleVDiff( xyz_UB, xyz_VB, xyz_TempB, xyr_Press, xyr_Temp, xyz_Height, xyz_Exner, xyr_Exner, xyz_DUDt, xyz_DVDt, xyz_DTempDt )

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, jmax, kmax    ! 鉛直層数.
                               ! Number of vertical level

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav, CpDry, GasRDry


    real(DP), intent(in ) :: xyz_UB     (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_VB     (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_TempB  (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyr_Press  (0:imax-1,1:jmax,0:kmax)
    real(DP), intent(in ) :: xyr_Temp   (0:imax-1,1:jmax,0:kmax)
    real(DP), intent(in ) :: xyz_Height (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyz_Exner  (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyr_Exner  (0:imax-1,1:jmax,0:kmax)
    real(DP), intent(out) :: xyz_DUDt   (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_DVDt   (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(out) :: xyz_DTempDt(0:imax-1,1:jmax,1:kmax)


    !
    ! local variables
    !
    real(dp)  :: xyr_tmdc (0:imax-1,1:jmax,0:kmax), xyr_thdc (0:imax-1,1:jmax,0:kmax)
    real(dp)  :: xyr_rho  (0:imax-1,1:jmax,0:kmax)
    real(dp)  :: xyr_tmflx(0:imax-1,1:jmax,0:kmax), xyr_tmfly(0:imax-1,1:jmax,0:kmax)
    real(dp)  :: xyr_thfl (0:imax-1,1:jmax,0:kmax)
    integer(4):: k


    do k = 0, 0
      xyr_tmdc(:,:,k) = 1.0d100
      xyr_thdc(:,:,k) = 1.0d100
    end do
    do k = 0+1, kmax-1
      xyr_tmdc(:,:,k) = VDiffCoefM
      xyr_thdc(:,:,k) = VDiffCoefH
    end do
    do k = kmax, kmax
      xyr_tmdc(:,:,k) = 1.0d100
      xyr_thdc(:,:,k) = 1.0d100
    end do

    do k = 0, 0
      xyr_rho(:,:,k) = 1.0d100
    end do
    do k = 0+1, kmax-1
      xyr_rho(:,:,k) = xyr_Press(:,:,k) / ( GasRDry * xyr_Temp(:,:,k) )
    end do
    do k = kmax, kmax
      xyr_rho(:,:,k) = 1.0d100
    end do

    do k = 0, 0
      xyr_tmflx(:,:,k) = 0.0d0
      xyr_tmfly(:,:,k) = 0.0d0
      xyr_thfl (:,:,k) = 0.0d0
    end do
    do k = 0+1, kmax-1
      xyr_tmflx(:,:,k) = - xyr_rho(:,:,k) * xyr_tmdc(:,:,k) * ( xyz_UB    (:,:,k+1) - xyz_UB    (:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
      xyr_tmfly(:,:,k) = - xyr_rho(:,:,k) * xyr_tmdc(:,:,k) * ( xyz_VB    (:,:,k+1) - xyz_VB    (:,:,k) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k) )
      xyr_thfl (:,:,k) = - CpDry * xyr_rho(:,:,k) * xyr_thdc(:,:,k) * xyr_Exner(:,:,k) * ( xyz_TempB(:,:,k+1) / xyz_Exner(:,:,k+1) - xyz_TempB(:,:,k  ) / xyz_Exner(:,:,k  ) ) / ( xyz_Height(:,:,k+1) - xyz_Height(:,:,k  ) )
    end do
    do k = kmax, kmax
      xyr_tmflx(:,:,k) = 0.0d0
      xyr_tmfly(:,:,k) = 0.0d0
      xyr_thfl (:,:,k) = 0.0d0
    end do


    do k = 1, kmax
      xyz_DUDt(:,:,k) = + Grav * ( xyr_tmflx(:,:,k) - xyr_tmflx(:,:,k-1) ) / ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) )
      xyz_DVDt(:,:,k) = + Grav * ( xyr_tmfly(:,:,k) - xyr_tmfly(:,:,k-1) ) / ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) )
      xyz_DTempDt(:,:,k) = + Grav / CpDry * ( xyr_thfl (:,:,k) - xyr_thfl (:,:,k-1) ) / ( xyr_Press(:,:,k) - xyr_Press(:,:,k-1) )
    end do


  end subroutine VenusSimpleVDiff