!= Yamamoto and Takahashi (2003) に従った簡単金星計算のための強制
!
!= forcing for simple Venus calculation following Yamamoto and Takahashi (2003)
!
! Authors::   Yoshiyuki O. Takahashi (code by Shin-ichi Takehiro is included)
! Version::   $Id: yt2003_forcing.f90,v 1.7 2025/09/17 22:00:00 takepiro Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2025. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module yt2003_forcing
  !
  != Yamamoto and Takahashi (2003) に従った簡単金星計算のための強制
  !
  != forcing for simple Venus calculation following Yamamoto and Takahashi (2003)

  !== Procedures List
  !
!!$  ! Hs94Forcing   :: 強制と散逸の計算
!!$  ! Hs94Finalize  :: 終了処理 (モジュール内部の変数の割り付け解除)
!!$  ! ------------  :: ------------
!!$  ! Hs94Forcing   :: Calculate forcing and dissipation
!!$  ! Hs94Finalize  :: Termination (deallocate variables in this module)
  !
  !--
  !== NAMELIST
  !
  ! NAMELIST#venus_simple_forcing_nml
  !++
  !== References
  !
  ! * Yamamoto, M, and M. Takahashi, 2003:
  !   The Fully Developed Superrotation Simulated by a General Circulation Model 
  !   of a Venus-like Atmosphere, 
  !   <i>J. Atmos. Sci.</i>, <b>60</b>, 561--574.
  ! * Hou, A. Y., and B. F. Farrell, 1987:
  !   Superrotation Induced by Critical-Level Absorption of Gravity Waves on Venus:
  !   An Assessment, 
  !   <i>J. Atmos. Sci.</i>, <b>44</b>, 1049--1061.
  !

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

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision.
    &                 STRING     ! 文字列.       Strings.

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

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  real(DP), parameter :: DayEarth = 86400.0d0

  real(DP), save      :: SurfFrictionTimeConstInEarthDay
  logical , save      :: FlagConstNCC
  real(DP), save      :: ConstNCCInEarthDay


  real(DP), save :: a_YT2003Temp      (0:25)
  real(DP), save :: a_YT2003HeightForT(0:25)
  real(DP), save :: a_YT2003Q         (0:20)
  real(DP), save :: a_YT2003HeightForQ(0:20)

  ! 公開手続き
  ! Public procedure
  !
  public :: YT2003Forcing
  public :: YT2003ForcingInit

  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !
  logical, save :: venus_simple_forcing_inited = .false.
                              ! 初期設定フラグ.
                              ! Initialization flag

  character(*), parameter:: module_name = 'venus_simple_forcing_1994'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: yt2003_forcing.f90,v 1.7 2025/09/17 22:00:00 takepiro Exp $'
                              ! モジュールのバージョン
                              ! Module version

  !--------------------------------------------------------------------------------------

contains

  subroutine YT2003Forcing(                                &
    & xy_SurfHeight,                                       & ! (in )
    & xyz_UB, xyz_VB, xyz_TempB, xyz_VirTemp, xyr_VirTemp, & ! (in )
    & xy_PsB, xyz_Press, xyr_Press, xyr_Temp,              & ! (in )
    & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner,        & ! (in )
    & xyz_DUDt, xyz_DVDt, xyz_DTempDt                      & ! (out)
    & )

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数.
                               ! Number of vertical level

    ! 組成に関わる配列の設定
    ! Settings of array for atmospheric composition
    !
    use composition, only: &
      &                    ncmax
                              ! 成分の数
                              ! Number of composition

    ! 鉛直拡散フラックス
    ! Vertical diffusion flux
    !
    use vdiffusion_my, only: VDiffusion, VDiffusionExpTendency

    real(DP), intent(in ) :: xy_SurfHeight(0:imax-1,1:jmax)
    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_VirTemp  (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyr_VirTemp  (0:imax-1,1:jmax,0:kmax)
    real(DP), intent(in ) :: xyz_Height   (0:imax-1,1:jmax,1:kmax)
    real(DP), intent(in ) :: xyr_Height   (0:imax-1,1:jmax,0: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_DTempDtRadL   (0:imax-1,1:jmax,1:kmax)
    real(DP) :: xyz_DTempDtRadS   (0:imax-1,1:jmax,1:kmax)

    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_DTempDtSFCFric(0:imax-1,1:jmax,1:kmax)

    real(DP) :: xyzf_QMix         (0:imax-1,1:jmax,1:kmax,1:ncmax)

    real(DP) :: xyr_MomFluxX      (0:imax-1,1:jmax,0:kmax)
    real(DP) :: xyr_MomFluxY      (0:imax-1,1:jmax,0:kmax)
    real(DP) :: xyr_HeatFlux      (0:imax-1,1:jmax,0:kmax)
    real(DP) :: xyrf_QMixFlux     (0:imax-1,1:jmax,0:kmax,1:ncmax)
    real(DP) :: xyr_VelTransCoef  (0:imax-1,1:jmax,0:kmax)
    real(DP) :: xyr_TempTransCoef (0:imax-1,1:jmax,0:kmax)
    real(DP) :: xyr_QMixTransCoef (0:imax-1,1:jmax,0: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)


    ! 初期化確認
    ! Initialization check
    !
    if ( .not. venus_simple_forcing_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    call YT2003RadForcing(           &
      & xy_PsB, xyz_Press, xyz_TempB, xyz_Height, & ! (in )
      & xyz_DTempDtRadL, xyz_DTempDtRadS          & ! (out)
      & )


    call YT2003SurfFriction( &
      & xyz_UB, xyz_VB, xyz_TempB,       & ! (in )
      & xyz_DUDtSFCFric, xyz_DVDtSFCFric, xyz_DTempDtSFCFric  & ! (out)
      & )


    ! This is set temporarily
    !
    xyzf_QMix = 0.0_DP

    call VDiffusion(                                               &
      & xyz_UB,     xyz_VB,     xyzf_QMix,                         & ! (in)
      & xyz_TempB, xyr_Temp, xyz_VirTemp, xyr_VirTemp, xyr_Press,  & ! (in)
      & xy_SurfHeight,                                             & ! (in)
      & xyz_Height, xyr_Height, xyz_Exner, xyr_Exner,              & ! (in)
      & xyr_MomFluxX,  xyr_MomFluxY,  xyr_HeatFlux, xyrf_QMixFlux, & ! (out)
      & xyr_VelTransCoef, xyr_TempTransCoef,                       & ! (out)
      & xyr_QMixTransCoef                                          & ! (out)
      & )

    call VDiffusionExpTendency(                            &
      & xyr_Press,                                               & ! (in )
      & xyr_MomFluxX, xyr_MomFluxY, xyr_HeatFlux, xyrf_QMixFlux, & ! (in ) optional
      & xyz_DUDtVDiff, xyz_DVDtVDiff, xyz_DTempDtVDiff  & ! (out) optional
      & )


    xyz_DUDt    =   xyz_DUDtVDiff    + xyz_DUDtSFCFric
    xyz_DVDt    =   xyz_DVDtVDiff    + xyz_DVDtSFCFric

    xyz_DTempDt =   xyz_DTempDtVDiff + xyz_DTempDtSFCFric &
      &           + xyz_DTempDtRadL  + xyz_DTempDtRadS


  end subroutine YT2003Forcing

  !--------------------------------------------------------------------------------------

  subroutine YT2003RadForcing(                &
    & xy_Ps, xyz_Press, xyz_Temp, xyz_Height, & ! (in )
    & xyz_DTempDtRadL, xyz_DTempDtRadS        & ! (out)
    & )

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $.
      & TimesetClockStart, TimesetClockStop

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数.
                               ! Number of vertical level

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, &               ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure
      & GasRDry
                              ! $ R $ [J kg-1 K-1].
                              ! 乾燥大気の気体定数.
                              ! Gas constant of air

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & y_Lat, &              ! $ \varphi $ [rad.] . 緯度. Latitude
      & 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  :: j, k
    integer  :: jp, jn


    y_CosLat = cos( y_Lat )

    call YT2003NCTempEq( &
      & xyz_Height, &
      & xyz_TempEq  &
      & )

    call YT2003NCCoef( &
      & xy_Ps, xyz_Press, &
      & xyz_NCC &
      & )

    xyz_DTempDtRadL = - xyz_NCC * ( xyz_Temp - xyz_TempEq )

    !
    !  add global mean cooling rate
    !
    xyz_DTempDtRadL = xyz_DTempDtRadL - xyz_YT2003Q0( xyz_Height ) / DayEarth


    call YT2003DTempDtRadS(   &
      & y_CosLat, xyz_Height,            & ! (in)
      & xyz_DTempDtRadS                  & ! (out)
      & )


    !
    ! code for debug
    !
!!$    do k = 1, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$           xyz_DTempDtRadS(i,j,k) = xyz_DTempDtRadS(i,j,k) &
!!$             & - Q0YT2003( xyz_Height(i,j,k) ) / DayEarth
!!$        end do
!!$      end do
!!$    end do
!!$
!!$    i = 0
!!$    do k = 1, kmax
!!$      do j = 1, jmax
!!$        write( 60, * ) j, xyz_Press(i,j,k), xyz_Height(i,j,k), xyz_DTempDtRadS(i,j,k) * DayEarth
!!$      end do
!!$      write( 60, * )
!!$    end do
!!$    call flush( 60 )
!!$
!!$    i = 0
!!$    j = jmax/2+1
!!$    do k = 1, kmax
!!$      write( 61, * ) k, xyz_Height(i,j,k), xyz_Press(i,j,k), &
!!$        & 1.0d0 / xyz_NCC(i,j,k) / DayEarth, xyz_TempEq(i,j,k), xyz_DTempDtRadS(i,j,k) * DayEarth
!!$    end do
!!$    call flush( 61 )
!!$    stop


    !
    ! In the following, output variables are calculated. 
    ! Variables calculated below are used only for output. 
    !

    xyz_EquivTempEq = xyz_TempEq                                      &
      & + ( xyz_DTempDtRadS - xyz_YT2003Q0( xyz_Height ) / DayEarth ) &
      &   / ( xyz_NCC + 1.0e-100_DP )


    ! 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.0_DP             &
      & - 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

    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 YT2003RadForcing

  !--------------------------------------------------------------------------------------

  subroutine YT2003NCTempEq( &
    & xyz_Height, & ! (in)
    & xyz_TempEq  & ! (out)
    & )

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: DP, &      ! 倍精度実数型. Double precision.
      &                 STRING     ! 文字列.       Strings.


    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数.
                               ! Number of vertical 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
    !

    xyz_TempEq = xyz_YT2003TempEq( xyz_Height )


  end subroutine YT2003NCTempEq

  !--------------------------------------------------------------------------------------

  subroutine YT2003NCCoef( &
    & xy_Ps, xyz_Press,         & ! (in)
    & xyz_NCC                   & ! (out)
    & )

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: DP, &      ! 倍精度実数型. Double precision.
      &                 STRING     ! 文字列.       Strings.


    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                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  :: i, j, k


    if ( FlagConstNCC ) then

      xyz_NCC = 1.0_DP / ( ConstNCCInEarthDay * DayEarth )

    else

      ! Thermal damping coefficient 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.0_DP ) then
              xyz_alp1(i,j,k) =  0.0_DP
              xyz_alp2(i,j,k) =  0.9_DP
              xyz_alp3(i,j,k) =  0.0_DP
            else if( -xyz_lnPRat(i,j,k) .le. 7.0_DP ) then
              xyz_alp1(i,j,k) = -4.5_DP
              xyz_alp2(i,j,k) =  2.0_DP
              xyz_alp3(i,j,k) =  5.0_DP
            else
              xyz_alp1(i,j,k) = -8.5_DP
              xyz_alp2(i,j,k) =  0.5_DP
              xyz_alp3(i,j,k) =  7.0_DP
            end if
          end do
        end do
      end do
      NCTimeConst0 = 1.32e9_DP
      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.0_DP / NCTimeConst
          end do
        end do
      end do

    end if


  end subroutine YT2003NCCoef

  !--------------------------------------------------------------------------------------

  subroutine YT2003DTempDtRadS(   &
    & y_CosLat, xyz_Height,            & ! (in)
    & xyz_DTempDtRadS                  & ! (out)
    & )

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: DP, &      ! 倍精度実数型. Double precision.
      &                 STRING     ! 文字列.       Strings.


    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数.
                               ! Number of vertical level

    use axesset, only: &
      & y_Lat_weight             ! $ \varphi $ [rad.] . 緯度重み. Latitude

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI
                              ! $ \pi $ .
                              ! 円周率.  Circular constant

    real(DP), intent(in ) :: y_CosLat       (1:jmax)
    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) :: GM
    integer  :: j
    integer  :: k


!!$    do k = 1, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$           xyz_DTempDtRadS(i,j,k) = Q0YT2003( xyz_Height(i,j,k) ) &
!!$             & * y_CosLat(j) * 4.0_DP / PI                        &
!!$             & / DayEarth
!!$        end do
!!$      end do
!!$    end do


    GM = sum( y_CosLat**(7.0_DP / 5.0_DP ) * y_Lat_weight ) / sum( y_Lat_Weight )

    xyz_DTempDtRadS = xyz_YT2003Q0( xyz_Height )
    do k = 1, kmax
      do j = 1, jmax
        xyz_DTempDtRadS(:,j,k) = xyz_DTempDtRadS(:,j,k)      &
          & * y_CosLat(j)**(7.0_DP/5.0_DP) / GM                    &
          & / DayEarth
      end do
    end do


  end subroutine YT2003DTempDtRadS

  !--------------------------------------------------------------------------------------

  subroutine YT2003SurfFriction(  &
    & xyz_UB, xyz_VB, xyz_TempB,       & ! (in )
    & xyz_DUDt, xyz_DVDt, xyz_DTempDt  & ! (out)
    & )

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

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数.
                               ! Number of vertical level

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, &               ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure


    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(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) :: SurfFrictionTimeConst
    real(DP) :: SurfTemp


    SurfFrictionTimeConst = SurfFrictionTimeConstInEarthDay * DayEarth

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

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


    xyz_DTempDt(:,:,2:kmax) = 0.0d0

    SurfTemp =                                                &
      &   ( a_YT2003Temp      (25) - A_YT2003Temp      (24) ) &
      & / ( a_YT2003HeightForT(25) - a_YT2003HeightForT(24) ) &
      & * ( 0.0_DP - a_YT2003HeightForT(24) )                 &
      & + A_YT2003Temp(24)
    xyz_DTempDt(:,:,1) = - ( xyz_TempB(:,:,1) - SurfTemp ) / SurfFrictionTimeConst


  end subroutine YT2003SurfFriction

  !--------------------------------------------------------------------------------------

  function xyz_YT2003TempEq( xyz_h_in )

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: DP         ! 倍精度実数型. Double precision.

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数.
                               ! Number of vertical level

    real(DP), intent(IN) :: xyz_h_in        (0:imax-1, 1:jmax, 1:kmax)     ! 高度(m)
    real(DP)             :: xyz_YT2003TempEq(0:imax-1, 1:jmax, 1:kmax)     ! 温度(K)

    !
    ! local variables
    !
    real(DP) :: x
    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: kk


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

          if ( xyz_h_in(i,j,k) > a_YT2003HeightForT(0) ) then
            xyz_YT2003TempEq(i,j,k) = a_YT2003Temp(0)
          else if ( xyz_h_in(i,j,k) < a_YT2003HeightForT(25) ) then
            xyz_YT2003TempEq(i,j,k) =                                 &
              &   ( a_YT2003Temp      (24) - a_YT2003Temp      (25) ) &
              & / ( a_YT2003HeightForT(24) - a_YT2003HeightForT(25) ) &
              & * ( xyz_h_in(i,j,k)        - a_YT2003HeightForT(25) ) &
              & + a_YT2003Temp(25)
          else
            do kk = 1, 25
              if ( ( xyz_h_in(i,j,k) < a_YT2003HeightForT(kk-1) ) .and. &
                &  ( xyz_h_in(i,j,k) > a_YT2003HeightForT(kk  ) ) ) then
                x =   ( a_YT2003HeightForT(kk-1) - xyz_h_in(i,j,k)        ) &
                  & / ( a_YT2003HeightForT(kk-1) - a_YT2003HeightForT(kk) )
                xyz_YT2003TempEq(i,j,k) = &
                  & ( 1 - x ) * a_YT2003Temp(kk-1) + x * a_YT2003Temp(kk)
              endif
            end do
          end if

        end do
      end do
    end do


  end function xyz_YT2003TempEq

  !--------------------------------------------------------------------------------------

  function xyz_YT2003Q0( xyz_h_in )

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: DP         ! 倍精度実数型. Double precision.

    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数.
                               ! Number of vertical level

    real(DP), intent(IN) :: xyz_h_in    (0:imax-1, 1:jmax, 1:kmax)         ! 高度(m)
    real(DP)             :: xyz_YT2003Q0(0:imax-1, 1:jmax, 1:kmax)         ! 温度(K)

    !
    ! local variables
    !
    real(DP) :: x
    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: kk


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

          if ( xyz_h_in(i,j,k) > a_YT2003HeightForQ(0) ) then
            xyz_YT2003Q0(i,j,k) = a_YT2003Q(0)
          else if ( xyz_h_in(i,j,k) <= a_YT2003HeightForQ(20) ) then
            xyz_YT2003Q0(i,j,k) = a_YT2003Q(20)
          else
            do kk = 1, 20
              if ( ( xyz_h_in(i,j,k) <= a_YT2003HeightForQ(kk-1) ) .and. &
                &  ( xyz_h_in(i,j,k) >  a_YT2003HeightForQ(kk  ) ) ) then
                x =   ( a_YT2003HeightForQ(kk-1) - xyz_h_in(i,j,k)        ) &
                  & / ( a_YT2003HeightForQ(kk-1) - a_YT2003HeightForQ(kk) )
                xyz_YT2003Q0(i,j,k) = ( 1 - x ) * a_YT2003Q(kk-1) + x * a_YT2003Q(kk)
              endif
            end do
          end if

        end do
      end do
    end do


  end function xyz_YT2003Q0

  !--------------------------------------------------------------------------------------

  subroutine YT2003ForcingInit
    !
    ! 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
    !

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


    ! 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

    ! 鉛直拡散フラックス
    ! Vertical diffusion flux
    !
    use vdiffusion_my, only: VDiffusionInit

    ! メッセージ制御
    ! Message control
    !
    use mpi_messagecntl, only : DoesOutputMPIMessage

    ! 宣言文 ; 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, &
      & FlagConstNCC,                    &
      & ConstNCCInEarthDay
          !
          ! デフォルト値については初期化手続 "venus_simple_forcing#YT2003ForcingInit"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "venus_simple_forcing#YT2003ForcingInit" for the default values.
          !

    ! 実行文 ; Executable statement
    !

    if ( venus_simple_forcing_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    SurfFrictionTimeConstInEarthDay = 30.0_DP
    FlagConstNCC                    = .false.
    ConstNCCInEarthDay              = 30.0_DP


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

      rewind( unit_nml )
      read( unit_nml,                     &  ! (in)
        & nml = venus_simple_forcing_nml, &  ! (out)
        & iostat = iostat_nml             &  ! (out)
        & )
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
      if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) 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' )


    a_YT2003Temp( 0)=178.76712328767124_DP ; a_YT2003HeightForT(0)=94.54545454545455e3_DP
    a_YT2003Temp( 1)=180.82191780821918_DP ; a_YT2003HeightForT(1)=89.35064935064935e3_DP
    a_YT2003Temp( 2)=184.93150684931507_DP ; a_YT2003HeightForT(2)=84.5021645021645e3_DP
    a_YT2003Temp( 3)=193.15068493150685_DP ; a_YT2003HeightForT(3)=80.0e3_DP
    a_YT2003Temp( 4)=203.4246575342466_DP  ; a_YT2003HeightForT(4)=75.84415584415585e3_DP
    a_YT2003Temp( 5)=217.8082191780822_DP  ; a_YT2003HeightForT(5)=71.34199134199135e3_DP
    a_YT2003Temp( 6)=236.3013698630137_DP  ; a_YT2003HeightForT(6)=67.18614718614718e3_DP
    a_YT2003Temp( 7)=252.73972602739724_DP ; a_YT2003HeightForT(7)=64.06926406926407e3_DP
    a_YT2003Temp( 8)=273.28767123287673_DP ; a_YT2003HeightForT(8)=60.60606060606061e3_DP
    a_YT2003Temp( 9)=295.8904109589041_DP  ; a_YT2003HeightForT(9)=56.45021645021645e3_DP
    a_YT2003Temp(10)=320.54794520547944_DP ; a_YT2003HeightForT(10)=52.98701298701299e3_DP
    a_YT2003Temp(11)=343.1506849315068_DP  ; a_YT2003HeightForT(11)=49.523809523809526e3_DP
    a_YT2003Temp(12)=365.75342465753425_DP ; a_YT2003HeightForT(12)=46.40692640692641e3_DP
    a_YT2003Temp(13)=392.4657534246575_DP  ; a_YT2003HeightForT(13)=42.5974025974026e3_DP
    a_YT2003Temp(14)=419.17808219178085_DP ; a_YT2003HeightForT(14)=38.78787878787879e3_DP
    a_YT2003Temp(15)=445.8904109589041_DP  ; a_YT2003HeightForT(15)=35.324675324675326e3_DP
    a_YT2003Temp(16)=472.6027397260274_DP  ; a_YT2003HeightForT(16)=31.861471861471863e3_DP
    a_YT2003Temp(17)=499.3150684931507_DP  ; a_YT2003HeightForT(17)=28.3982683982684e3_DP
    a_YT2003Temp(18)=528.0821917808219_DP  ; a_YT2003HeightForT(18)=24.935064935064936e3_DP
    a_YT2003Temp(19)=556.8493150684931_DP  ; a_YT2003HeightForT(19)=21.125541125541126e3_DP
    a_YT2003Temp(20)=587.6712328767123_DP  ; a_YT2003HeightForT(20)=17.316017316017316e3_DP
    a_YT2003Temp(21)=614.3835616438356_DP  ; a_YT2003HeightForT(21)=13.852813852813853e3_DP
    a_YT2003Temp(22)=645.2054794520548_DP  ; a_YT2003HeightForT(22)=10.043290043290042e3_DP
    a_YT2003Temp(23)=669.8630136986301_DP  ; a_YT2003HeightForT(23)=6.926406926406926e3_DP
    a_YT2003Temp(24)=698.6301369863014_DP  ; a_YT2003HeightForT(24)=3.463203463203463e3_DP
    a_YT2003Temp(25)=725.3424657534247_DP  ; a_YT2003HeightForT(25)=0.3463203463203463e3_DP


    a_YT2003Q( 0)=0.0_DP                ; a_YT2003HeightForQ( 0)=80.0e3_DP
    a_YT2003Q( 1)=0.1282051282051282_DP ; a_YT2003HeightForQ( 1)=73.33333333333333e3_DP
    a_YT2003Q( 2)=0.2564102564102564_DP ; a_YT2003HeightForQ( 2)=71.11111111111111e3_DP
    a_YT2003Q( 3)=0.7692307692307693_DP ; a_YT2003HeightForQ( 3)=68.33333333333333e3_DP
    a_YT2003Q( 4)=1.4102564102564104_DP ; a_YT2003HeightForQ( 4)=66.11111111111111e3_DP
    a_YT2003Q( 5)=2.051282051282051_DP  ; a_YT2003HeightForQ( 5)=64.44444444444444e3_DP
    a_YT2003Q( 6)=2.6923076923076925_DP ; a_YT2003HeightForQ( 6)=62.77777777777778e3_DP
    a_YT2003Q( 7)=3.3333333333333335_DP ; a_YT2003HeightForQ( 7)=61.666666666666664e3_DP
    a_YT2003Q( 8)=3.9743589743589745_DP ; a_YT2003HeightForQ( 8)=60.55555555555556e3_DP
    a_YT2003Q( 9)=4.743589743589744_DP  ; a_YT2003HeightForQ( 9)=58.333333333333336e3_DP
    a_YT2003Q(10)=5.128205128205129_DP  ; a_YT2003HeightForQ(10)=56.666666666666664e3_DP
    a_YT2003Q(11)=5.2_DP                ; a_YT2003HeightForQ(11)=55.0e3_DP
    a_YT2003Q(12)=4.871794871794871_DP  ; a_YT2003HeightForQ(12)=52.22222222222222e3_DP
    a_YT2003Q(13)=4.358974358974359_DP  ; a_YT2003HeightForQ(13)=50.0e3_DP
    a_YT2003Q(14)=3.58974358974359_DP   ; a_YT2003HeightForQ(14)=47.77777777777778e3_DP
    a_YT2003Q(15)=2.948717948717949_DP  ; a_YT2003HeightForQ(15)=46.111111111111114e3_DP
    a_YT2003Q(16)=2.3076923076923075_DP ; a_YT2003HeightForQ(16)=43.888888888888886e3_DP
    a_YT2003Q(17)=1.6666666666666667_DP ; a_YT2003HeightForQ(17)=41.666666666666664e3_DP
    a_YT2003Q(18)=1.1538461538461537_DP ; a_YT2003HeightForQ(18)=38.888888888888886e3_DP
    a_YT2003Q(19)=0.7692307692307693_DP ; a_YT2003HeightForQ(19)=37.22222222222222e3_DP
    a_YT2003Q(20)=0.52_DP               ; a_YT2003HeightForQ(20)=35.0e3_DP


    ! Initialization of modules used in this module
    !

    ! 鉛直拡散フラックス (Mellor and Yamada, 1974)
    ! Vertical diffusion flux (Mellor and Yamada, 1974)
    !
    call VDiffusionInit


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, '  SurfFrictionTimeConstInEarthDay = %f', d = (/ SurfFrictionTimeConstInEarthDay /) )
    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 YT2003ForcingInit

  !--------------------------------------------------------------------------------------

  !
  ! A subroutine below will be deleted (yot, 2010/10/29)
  !
  subroutine VenusSimpleNCTempEq_old( &
    & xyz_Height, & ! (in)
    & xyz_TempEq  & ! (out)
    & )

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: DP, &      ! 倍精度実数型. Double precision.
      &                 STRING     ! 文字列.       Strings.


    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                kmax    ! 鉛直層数.
                               ! Number of vertical level

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, &               ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & CpDry
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure


    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.0e3_DP
    z ( 2 ) =  10.0e3_DP
    z ( 3 ) =  25.0e3_DP
!!$    z ( 4 ) =  55.0e3_DP
    z ( 4 ) =  50.0e3_DP
    z ( 5 ) = 100.0e3_DP

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

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



    a ( 1 ) =   0.0e0_DP

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


    SurfTemp = 750.0_DP
    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.5_DP &
        &   * ( 1.0_DP + tanh( ( 0.0_DP      - z(l) ) / d(l) ) )
      xyz_TempEq = xyz_TempEq &
        & + ( a(l+1) - a(l) ) * 0.5_DP &
        &   * ( 1.0_DP + 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_old

  !--------------------------------------------------------------------------------------

  !
  ! A subroutine below will be deleted (yot, 2010/10/29)
  !

  subroutine VenusSimpleDTempDtRadS_old(   &
    & y_CosLat, xyz_Press, xyz_Height, & ! (in)
    & xyz_DTempDtRadS                  & ! (out)
    )

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav, &               ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration
      & CpDry, &
                              ! $ C_p $ [J kg-1 K-1].
                              ! 乾燥大気の定圧比熱.
                              ! Specific heat of air at constant pressure
      & GasRDry

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: DP, &      ! 倍精度実数型. Double precision.
      &                 STRING     ! 文字列.       Strings.


    ! 格子点設定
    ! Grid points settings
    !
    use gridset, only: imax, & ! 経度格子点数.
                               ! Number of grid points in longitude
      &                jmax, & ! 緯度格子点数.
                               ! Number of grid points in latitude
      &                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
    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.0_DP / Grav

    xyz_DTempDtRadS &
      & = 5.0_DP / DayEarth &
      & * exp( - ( ( - scaleheight * log( xyz_Press / 500.0e2_DP ) ) / ( 2.0_DP * scaleheight ) )**2 )

    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1
          if ( xyz_Press(i,j,k) > 500.0e2_DP ) then
            if ( xyz_DTempDtRadS(i,j,k) .lt. 0.5_DP / DayEarth ) then
              xyz_DTempDtRadS(i,j,k) = 0.5_DP / 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_old

  !--------------------------------------------------------------------------------------

end module yt2003_forcing