| Class | hs94forcing_mod |
| In: |
hs94forcing/hs94forcing.f90
|
Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!=begin
* Developers: Morikawa Yasuhiro * Version: $Id: hs94forcing.f90,v 1.9 2005/01/22 09:23:41 morikawa Exp $ * Tag Name: $Name: $ * Change History:
This module compute heating and dissipation for Held and Suarez(1994) benchmark integration of a dry GCM. Held and Suarez(1994) の乾燥大気 GCM ベンチマーク用の加熱と散逸を計算する。
A proposal for the intercomparison of the dynamical cores of atmospheric general circuation models. Bull. Am. Meteor. Soc., 75, 1825--1830.
に依存するようになっている。正しくテストしたいのならば、それらが Held and Suarez(1994) の実験設定と整合的かどうかチェックするように すべきである。
!=end
加熱と散逸の結果は xyz_VelLon_phy, xyz_VelLat_phy, xyz_Temp_phy として返る。
subroutine hs94forcing (
xyz_VelLon_b , xyz_VelLat_b , xyz_Temp_b , xy_Ps_b ,
xyz_VelLon_phy, xyz_VelLat_phy, xyz_Temp_phy )
!==== Dependency
use type_mod, only: INTKIND, STRING, TOKEN, REKIND, DBKIND
use constants_mod, only: RAir, Cp, SecPerDay
use grid_3d_mod , only: im, jm, km
use io_gt4_out_mod,only: io_gt4_out_Put
use dc_trace , only: DbgMessage, BeginSub, EndSub, DataDump
!=end
implicit none
!=begin
!==== Input
!
real(DBKIND), intent(in) :: &
& xyz_VelLon_b(:,:,:) , & ! 速度経度成分 (t-Δt)
& xyz_VelLat_b(:,:,:) , & ! 速度緯度成分 (t-Δt)
& xyz_Temp_b(:,:,:) , & ! 温度 (t-Δt)
& xy_Ps_b(:,:) ! 地表面気圧 (t-Δt)
!==== Output
!
real(DBKIND), intent(out) :: &
& xyz_VelLon_phy(:,:,:) , & ! 速度経度成分の加熱散逸効果
& xyz_VelLat_phy(:,:,:) , & ! 速度緯度成分の加熱散逸効果
& xyz_Temp_phy(:,:,:) ! 温度の加熱散逸効果
!=end
!----------------------------------------------------------------
! 全体で用いる変数
!----------------------------------------------------------------
real(DBKIND) :: SigmaB ! σ_b
!----------------------------------------------------------------
! xyz_VelLon_phy, xyz_VelLat_phy を求めるための変数
!----------------------------------------------------------------
real(DBKIND) :: kf ! k_f
!----------------------------------------------------------------
! xyz_Temp_phy を求めるための変数
!----------------------------------------------------------------
real(DBKIND) :: &
& Kappa , & ! κ = Ｒ/Ｃp
& Press0 , & ! Ｐ_0 (Pa)
& ka , & ! k_a
& ks , & ! k_s
& DelTempY , & ! (ΔT)_y
& DelPotTempZ ! (Δθ)_z
!----------------------------------------------------------------
! 汎用変数
!----------------------------------------------------------------
integer(INTKIND) :: i, j, k
character(STRING), parameter:: subname = "hs94forcing"
continue
!----------------------------------------------------------------
! Check Initialization
!----------------------------------------------------------------
call BeginSub(subname)
if (.not. hs94forcing_initialized) then
call EndSub( subname, 'Call hs94forcing_init before call %c', &
& c1=trim(subname) )
return
endif
!----------------------------------------------------------------
! 全体で用いる変数
!----------------------------------------------------------------
SigmaB = 0.7d0
!----------------------------------------------------------------
! xyz_VelLon_phy, xyz_VelLat_phy を求める
!----------------------------------------------------------------
kf = 1.0d0 / SecPerDay
xyz_kv = kf * max( 0.0d0, ( xyz_Sigma - SigmaB ) /( 1.0d0 - SigmaB ) )
xyz_VelLon_phy = - xyz_kv * xyz_VelLon_b
xyz_VelLat_phy = - xyz_kv * xyz_VelLat_b
call DbgMessage('kf=<%f>', d=(/kf/))
call io_gt4_out_Put('xyz_VelLon_phy', xyz_VelLon_phy)
call io_gt4_out_Put('xyz_VelLat_phy', xyz_VelLat_phy)
call io_gt4_out_Put('xyz_kv', xyz_kv)
!----------------------------------------------------------------
! xyz_Temp_phy を求める
!----------------------------------------------------------------
do k = 1, km
xyz_Press_b(:,:,k) = xyz_Sigma(:,:,k) * xy_Ps_b(:,:)
enddo
Kappa = RAir / Cp ! κ=Ｒ／Ｃｐ (気体定数/定圧比熱)
Press0 = 1000.0d0 * 1.0d2
ka = 1.0d0 / ( 40.0d0 * SecPerDay )
ks = 1.0d0 / ( 4.0d0 * SecPerDay )
DelTempY = 60.0d0
DelPotTempZ = 10.0d0
xyz_TempEQ = &
& max( 200.0d0, &
& ( 315.0d0 - DelTempY * xyz_SinLat**2 &
& - DelPotTempZ * log( xyz_Press_b / Press0 ) &
& * xyz_CosLat**2 &
& ) &
& * ( xyz_Press_b / Press0 )**Kappa &
& )
xyz_kt = ka + ( ks - ka ) &
& * max( 0.0d0, &
& ( xyz_Sigma - SigmaB ) / ( 1.0d0 - SigmaB ) &
& ) &
& * xyz_CosLat**4
xyz_Temp_phy = - xyz_kt * ( xyz_Temp_b - xyz_TempEQ )
call io_gt4_out_Put('xyz_Temp_phy', xyz_Temp_phy)
call io_gt4_out_Put('xyz_TempEQ', xyz_TempEQ)
call io_gt4_out_Put('xyz_PotTempEQ', &
& xyz_TempEQ * ( xyz_Sigma )**( -Kappa ) )
!!$ call io_gt4_out_Put('1.0/max(xyz_kv, 1.0/86400d2)/86400.0', &
!!$ & 1.0/max(xyz_kv, 1.0/86400d2)/86400.0 )
!!$ call io_gt4_out_Put('1.0/xyz_kt', 1.0/xyz_kt)
!!$ call DataDump('1.0/xyz_kt', 1.0/xyz_kt, strlen=60)
!!$
!!$ call DataDump('xyz_Press_b', xyz_Press_b, strlen=60)
!!$
!!$ call DataDump('xyz_TempEQ1', &
!!$ & ( 315.0d0 - DelTempY * xyz_SinLat**2 &
!!$ & - DelPotTempZ * log( xyz_Press_b / Press0 ) &
!!$ & * xyz_CosLat**2 &
!!$ & ) &
!!$ & * ( xyz_Press_b / Press0 )**Kappa &
!!$ & , strlen=60 )
!!$
!!$ call DataDump('xyz_TempEQ2', &
!!$ & - DelTempY * xyz_SinLat**2 &
!!$ & , strlen=60 )
!!$
!!$ call DataDump('xyz_TempEQ3', &
!!$ & ( - DelPotTempZ * log( xyz_Press_b / Press0 ) &
!!$ & * xyz_CosLat**2 &
!!$ & ) &
!!$ & , strlen=60 )
call EndSub(subname)
end subroutine
以降のサブルーチンで用いる変数の allocate を行なう。
本来はここで、実験設定と ((< constants_mod >)) の値との比較を 行ないたいが、まだ実装されていない。 (実数同士のそのままの比較は難しいので、有効数字数桁と指数部分を 比較するようにすべきだろう)。
subroutine hs94forcing_init (
x_Lon , y_Lat , z_Sigma )
!==== Dependency
use type_mod, only: INTKIND, STRING, TOKEN, REKIND, DBKIND
use constants_mod, only: constants_init, pi
use grid_3d_mod , only: grid_3d_init, im, jm, km
use axis_type_mod, only: AXISINFO
use io_gt4_out_mod,only: io_gt4_out_init, io_gt4_out_SetVars
use dc_string , only: LChar, StrHead
use dc_trace , only: DbgMessage, BeginSub, EndSub, DataDump
!=end
implicit none
!=begin
!==== Input
!
type(AXISINFO), intent(in) :: &
& x_Lon , & ! 経度座標
& y_Lat , & ! 緯度座標
& z_Sigma ! σレベル(整数)座標
!=end
integer(INTKIND) :: i, j, k
real(DBKIND) :: RadDegFact ! ラジアンと度数の変換係数
character(STRING), parameter:: subname = "hs94forcing_init"
continue
!----------------------------------------------------------------
! Check Initialization
!----------------------------------------------------------------
call BeginSub(subname)
if (hs94forcing_initialized) then
call EndSub( subname, '%c is already called.', c1=trim(subname) )
return
else
hs94forcing_initialized = .true.
endif
!----------------------------------------------------------------
! Version identifier
!----------------------------------------------------------------
call DbgMessage('%c :: %c', c1=trim(version), c2=trim(tagname))
!-----------------------------------------------------------------
! Initialize Dependent modules
!-----------------------------------------------------------------
call constants_init
call grid_3d_init
call io_gt4_out_init
!-------------------------------------------------------------------
! Setting Ouput Data by io_gt4_out_set_Vars
!-------------------------------------------------------------------
call io_gt4_out_SetVars('xyz_VelLon_phy')
call io_gt4_out_SetVars('xyz_VelLat_phy')
call io_gt4_out_SetVars('xyz_Temp_phy')
call io_gt4_out_SetVars('xyz_TempEQ')
call io_gt4_out_SetVars('xyz_PotTempEQ')
call io_gt4_out_SetVars('1.0/xyz_kt')
call io_gt4_out_SetVars('1.0/max(xyz_kv, 1.0/86400d2)/86400.0')
call io_gt4_out_SetVArs('xyz_kv')
!-----------------------------------------------------------------
! Allocate variables
!-----------------------------------------------------------------
allocate( &
& xyz_Sigma(im, jm, km) , & ! σ
& xyz_kv(im, jm, km) , & ! k_v
& xyz_Press_b(im, jm, km) , & ! 気圧
& xyz_SinLat(im, jm, km) , & ! sinφ
& xyz_CosLat(im, jm, km) , & ! cosφ
& xyz_TempEQ(im,jm,km) , & ! T_eq
& xyz_kt(im,jm,km) ) ! k_t
!----------------------------------------------------------------
! 全体で用いる変数
!----------------------------------------------------------------
do k = 1, km
xyz_Sigma(:,:,k) = z_Sigma%a_Dim(k)
enddo
!----------------------------------------------------------------
! xyz_Temp_phy を求めるための変数
!----------------------------------------------------------------
if ( StrHead( 'radians', trim(LChar(y_Lat%axisinfo%units)) ) .or.&
& StrHead( 'rad.', trim(LChar(y_Lat%axisinfo%units)) ) ) then
RadDegFact = 1.
else
RadDegFact = pi / 180.
end if
do j = 1, jm
xyz_SinLat(:,j,:) = sin( y_Lat%a_Dim(j) * RadDegFact )
xyz_CosLat(:,j,:) = cos( y_Lat%a_Dim(j) * RadDegFact )
enddo
call DataDump('xyz_SinLat', xyz_SinLat, strlen=70)
call DataDump('xyz_CosLat', xyz_CosLat, strlen=70)
call EndSub(subname)
end subroutine