!= 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.1 2010-11-14 12:36:19 yot Exp $
! Tag Name:: $Name: dcpam5-20110615 $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. 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,
! J. Atmos. Sci., 60, 561--574.
! * Hou, A. Y., and B. F. Farrell, 1987:
! Superrotation Induced by Critical-Level Absorption of Gravity Waves on Venus:
! An Assessment,
! J. Atmos. Sci., 44, 1049--1061.
!
! モジュール引用 ; USE statements
!
! 種別型パラメタ
! Kind type parameter
!
use dc_types, only: DP, & ! 倍精度実数型. Double precision.
& STRING ! 文字列. Strings.
! 宣言文 ; 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 variables
!
logical, save, public:: venus_simple_forcing_inited = .false.
! 初期設定フラグ.
! Initialization flag
character(*), parameter:: module_name = 'venus_simple_forcing_1994'
! モジュールの名称.
! Module name
character(*), parameter:: version = &
& '$Name: dcpam5-20110615 $' // &
& '$Id: yt2003_forcing.f90,v 1.1 2010-11-14 12:36:19 yot Exp $'
! モジュールのバージョン
! Module version
!--------------------------------------------------------------------------------------
contains
subroutine YT2003Forcing( &
& xy_SurfHeight, & ! (in )
& xyz_UB, xyz_VB, xyz_TempB, & ! (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
! 鉛直拡散フラックス (Mellor and Yamada, 1974, レベル 2)
! Vertical diffusion flux (Mellor and Yamada, 1974, Level 2)
!
use vdiffusion_my1974, 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_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
!
if ( .not. venus_simple_forcing_inited ) call YT2003ForcingInit
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, 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: &
& DelTime, & ! $ \Delta t $
& 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.0d-100 )
! 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
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.0d0 / ( 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.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
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 constants settings
!
use constants, 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/5.0) / 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
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(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) :: Temp (24:25)
real(DP) :: Height(24:25)
real(DP) :: SurfTemp
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
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.
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! 格子点設定
! 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.
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! 格子点設定
! 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
!
! メッセージ出力
! Message output
!
use dc_message, only: MessageNotify
! 物理定数設定
! Physical constants settings
!
use constants, only: &
& GasRDry, &
! $ R $ [J kg-1 K-1].
! 乾燥大気の気体定数.
! Gas constant of air
& CpDry
! $ C_p $ [J kg-1 K-1].
! 乾燥大気の定圧比熱.
! Specific heat of air at constant pressure
! 座標データ設定
! Axes data settings
!
use axesset, only: &
& y_Lat, & ! $ \varphi $ [rad.] . 緯度. Latitude
& 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, &
& 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.0d0
FlagConstNCC = .false.
ConstNCCInEarthDay = 30.0d0
! 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 ) 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 ; a_YT2003HeightForT(0)=94.54545454545455d3
a_YT2003Temp( 1)=180.82191780821918 ; a_YT2003HeightForT(1)=89.35064935064935d3
a_YT2003Temp( 2)=184.93150684931507 ; a_YT2003HeightForT(2)=84.5021645021645d3
a_YT2003Temp( 3)=193.15068493150685 ; a_YT2003HeightForT(3)=80.0d3
a_YT2003Temp( 4)=203.4246575342466 ; a_YT2003HeightForT(4)=75.84415584415585d3
a_YT2003Temp( 5)=217.8082191780822 ; a_YT2003HeightForT(5)=71.34199134199135d3
a_YT2003Temp( 6)=236.3013698630137 ; a_YT2003HeightForT(6)=67.18614718614718d3
a_YT2003Temp( 7)=252.73972602739724 ; a_YT2003HeightForT(7)=64.06926406926407d3
a_YT2003Temp( 8)=273.28767123287673 ; a_YT2003HeightForT(8)=60.60606060606061d3
a_YT2003Temp( 9)=295.8904109589041 ; a_YT2003HeightForT(9)=56.45021645021645d3
a_YT2003Temp(10)=320.54794520547944 ; a_YT2003HeightForT(10)=52.98701298701299d3
a_YT2003Temp(11)=343.1506849315068 ; a_YT2003HeightForT(11)=49.523809523809526d3
a_YT2003Temp(12)=365.75342465753425 ; a_YT2003HeightForT(12)=46.40692640692641d3
a_YT2003Temp(13)=392.4657534246575 ; a_YT2003HeightForT(13)=42.5974025974026d3
a_YT2003Temp(14)=419.17808219178085 ; a_YT2003HeightForT(14)=38.78787878787879d3
a_YT2003Temp(15)=445.8904109589041 ; a_YT2003HeightForT(15)=35.324675324675326d3
a_YT2003Temp(16)=472.6027397260274 ; a_YT2003HeightForT(16)=31.861471861471863d3
a_YT2003Temp(17)=499.3150684931507 ; a_YT2003HeightForT(17)=28.3982683982684d3
a_YT2003Temp(18)=528.0821917808219 ; a_YT2003HeightForT(18)=24.935064935064936d3
a_YT2003Temp(19)=556.8493150684931 ; a_YT2003HeightForT(19)=21.125541125541126d3
a_YT2003Temp(20)=587.6712328767123 ; a_YT2003HeightForT(20)=17.316017316017316d3
a_YT2003Temp(21)=614.3835616438356 ; a_YT2003HeightForT(21)=13.852813852813853d3
a_YT2003Temp(22)=645.2054794520548 ; a_YT2003HeightForT(22)=10.043290043290042d3
a_YT2003Temp(23)=669.8630136986301 ; a_YT2003HeightForT(23)=6.926406926406926d3
a_YT2003Temp(24)=698.6301369863014 ; a_YT2003HeightForT(24)=3.463203463203463d3
a_YT2003Temp(25)=725.3424657534247 ; a_YT2003HeightForT(25)=0.3463203463203463d3
a_YT2003Q( 0)=0.0 ; a_YT2003HeightForQ( 0)=80.0d3
a_YT2003Q( 1)=0.1282051282051282 ; a_YT2003HeightForQ( 1)=73.33333333333333d3
a_YT2003Q( 2)=0.2564102564102564 ; a_YT2003HeightForQ( 2)=71.11111111111111d3
a_YT2003Q( 3)=0.7692307692307693 ; a_YT2003HeightForQ( 3)=68.33333333333333d3
a_YT2003Q( 4)=1.4102564102564104 ; a_YT2003HeightForQ( 4)=66.11111111111111d3
a_YT2003Q( 5)=2.051282051282051 ; a_YT2003HeightForQ( 5)=64.44444444444444d3
a_YT2003Q( 6)=2.6923076923076925 ; a_YT2003HeightForQ( 6)=62.77777777777778d3
a_YT2003Q( 7)=3.3333333333333335 ; a_YT2003HeightForQ( 7)=61.666666666666664d3
a_YT2003Q( 8)=3.9743589743589745 ; a_YT2003HeightForQ( 8)=60.55555555555556d3
a_YT2003Q( 9)=4.743589743589744 ; a_YT2003HeightForQ( 9)=58.333333333333336d3
a_YT2003Q(10)=5.128205128205129 ; a_YT2003HeightForQ(10)=56.666666666666664d3
a_YT2003Q(11)=5.2 ; a_YT2003HeightForQ(11)=55.0d3
a_YT2003Q(12)=4.871794871794871 ; a_YT2003HeightForQ(12)=52.22222222222222d3
a_YT2003Q(13)=4.358974358974359 ; a_YT2003HeightForQ(13)=50.0d3
a_YT2003Q(14)=3.58974358974359 ; a_YT2003HeightForQ(14)=47.77777777777778d3
a_YT2003Q(15)=2.948717948717949 ; a_YT2003HeightForQ(15)=46.111111111111114d3
a_YT2003Q(16)=2.3076923076923075 ; a_YT2003HeightForQ(16)=43.888888888888886d3
a_YT2003Q(17)=1.6666666666666667 ; a_YT2003HeightForQ(17)=41.666666666666664d3
a_YT2003Q(18)=1.1538461538461537 ; a_YT2003HeightForQ(18)=38.888888888888886d3
a_YT2003Q(19)=0.7692307692307693 ; a_YT2003HeightForQ(19)=37.22222222222222d3
a_YT2003Q(20)=0.52 ; a_YT2003HeightForQ(20)=35.0d3
! 印字 ; 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
use axesset, only: &
& y_Lat, & ! $ \varphi $ [rad.] . 緯度. Latitude
& 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_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
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_old
!--------------------------------------------------------------------------------------
end module yt2003_forcing