| Class | dynamics_mod |
| In: |
dynamics/dynamics.f90
|
| Developers : | Morikawa Yasuhiro |
| Version : | $Id: dynamics.f90,v 1.20 2005/01/22 09:30:59 morikawa Exp $ |
| Tag Name : | $Name: $ |
Change History ::
Calculate Dynamical Core. 力学コア部分を演算するモジュール。演算している方程式系、 および離散化の手法は以下の通り。
特に無し。
特に無し
現在、力学コアの全てがこのモジュール内にある。本来ならば、 解く方程式や項の種類によってモジュール化がなされるべきである。
診断的に得られる量の演算を行なう。 現在は渦度発散から速度成分 (経度方向、緯度方向) の演算を行なうのみである。
subroutine dynamics_diagnostic (
x_Lon , y_Lat , z_Sigma , r_Sigma ,
xyz_VelLon_a, xyz_VelLat_a, xyz_Vor_a, xyz_Div_a ,
xyz_Temp_a , xyz_QVap_a , xy_Ps_a )
!==== Dependency
use type_mod, only: STRING, REKIND, DBKIND, INTKIND
use grid_3d_mod, only: im, jm, km
use grid_wavenumber_mod, only: nm
use constants_mod, only: R0, Grav
use spml_mod, only: wa_xya, xya_GradLon_wa, xya_GradLat_wa, &
& wa_LaplaInv_wa, IntLonLat_xy
use io_gt4_out_mod,only: io_gt4_out_Put
use dc_trace, only: DbgMessage, BeginSub, EndSub, DataDump
use dc_string, only: toChar
implicit none
!==== Input
!
real(DBKIND), intent(in) :: &
& x_Lon(:) , & ! intent(in): 経度座標
& y_Lat(:) , & ! intent(in): 緯度座標
& z_Sigma(:) , & ! intent(in): σレベル(整数)座標
& r_Sigma(:) , & ! intent(in): σレベル(半整数)座標
&
& xyz_Vor_a(:,:,:) , & ! intent(in): 渦度 (t+Δt)
& xyz_Div_a(:,:,:) , & ! intent(in): 発散 (t+Δt)
& xyz_Temp_a(:,:,:) , & ! intent(in): 温度 (t+Δt)
& xyz_QVap_a(:,:,:) , & ! intent(in): 比湿 (t+Δt)
& xy_Ps_a(:,:) ! intent(in): 地表面気圧 (t+Δt)
!==== In/Output
!
real(DBKIND), intent(out) :: &
& xyz_VelLon_a(:,:,:) , & ! intent(out): 速度経度成分 (t+Δt)
& xyz_VelLat_a(:,:,:) ! intent(out): 速度緯度成分 (t+Δt)
real(DBKIND) :: TotalMass ! 全質量
!----- 作業用内部変数 -----
! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
integer(INTKIND) :: i, j, k, l
character(STRING), parameter:: subname = "dynamics_diagnostic"
!----------------------------------------------------------------
! Check Initialization
!----------------------------------------------------------------
call BeginSub(subname)
if (.not. dynamics_initialized) then
call EndSub( subname, 'Call dynamics_init before call %c', &
& c1=trim(subname) )
return
endif
!-------------------------------------------------------------------
! ｕ と ｖ を出力
!-------------------------------------------------------------------
wz_Psi_a(:,:) = wa_LaplaInv_wa( wa_xya( xyz_Vor_a(:,:,:) ) ) * R0**2
wz_Chi_a(:,:) = wa_LaplaInv_wa( wa_xya( xyz_Div_a(:,:,:) ) ) * R0**2
xyz_VelLon_a(:,:,:) = ( xya_GradLon_wa( wz_Chi_a(:,:) ) &
& - xya_GradLat_wa( wz_Psi_a(:,:) ) ) / R0
xyz_VelLat_a(:,:,:) = ( xya_GradLon_wa( wz_Psi_a(:,:) ) &
& + xya_GradLat_wa( wz_Chi_a(:,:) ) ) / R0
!!$ call DataDump('xyz_VelLon_a', xyz_VelLon_a(:,:,:), strlen=80)
!!$ call DataDump('xyz_VelLat_a', xyz_VelLat_a(:,:,:), strlen=80)
!-------------------------------------------------------------------
! xy_Ps から全質量を求める。
!-------------------------------------------------------------------
TotalMass = IntLonLat_xy( xy_Ps_a / Grav )
call io_gt4_out_Put('TotalMass', TotalMass)
call EndSub(subname)
end subroutine
t-Δt の値から水平拡散項を求め、それを t+Δt の値に加える。
subroutine dynamics_diffusion (
xyz_Vor_b , xyz_Div_b , xyz_Temp_b , xyz_QVap_b ,
xyz_Vor_a , xyz_Div_a , xyz_Temp_a , xyz_QVap_a )
!==== Dependency
use type_mod, only: STRING, REKIND, DBKIND, INTKIND
use time_mod, only: DelTime, CurrentTime
use grid_3d_mod, only: km
use grid_wavenumber_mod, only: nm
use spml_mod, only: wa_xya, xya_wa, l_nm
use io_gt4_out_mod,only: io_gt4_out_Put
use dc_trace, only: DbgMessage, BeginSub, EndSub, DataDump
implicit none
!==== Input
!
real(DBKIND), intent(in) :: &
& xyz_Vor_b(:,:,:) , & ! 渦度 (t-Δt)
& xyz_Div_b(:,:,:) , & ! 発散 (t-Δt)
& xyz_Temp_b(:,:,:) , & ! 温度 (t-Δt)
& xyz_QVap_b(:,:,:) ! 比湿 (t-Δt)
!==== In/Out
!
real(DBKIND), intent(inout) :: &
& xyz_Vor_a(:,:,:) , & ! 渦度 (t+Δt)
& xyz_Div_a(:,:,:) , & ! 発散 (t+Δt)
& xyz_Temp_a(:,:,:) , & ! 温度 (t+Δt)
& xyz_QVap_a(:,:,:) ! 比湿 (t+Δt)
real(DBKIND) :: &
& wz_Vor_a((nm+1)**2, km) , & ! 渦度 (t+Δt)
& wz_Div_a((nm+1)**2, km) , & ! 発散 (t+Δt)
& wz_Temp_a((nm+1)**2, km) , & ! 温度 (t+Δt)
& wz_QVap_a((nm+1)**2, km) ! 比湿 (t+Δt)
integer(INTKIND) :: i
character(STRING), parameter:: subname = "dynamics_diffusion"
!----------------------------------------------------------------
! Check Initialization
!----------------------------------------------------------------
call BeginSub(subname)
if (.not. dynamics_initialized) then
call EndSub( subname, 'Call dynamics_init before call %c', &
& c1=trim(subname) )
return
endif
!-------------------------------------------------------------------
! それぞれの量の拡散項を計算
!-------------------------------------------------------------------
xyz_Vor_a = &
& xya_wa( &
& wa_xya(xyz_Vor_a) &
& + 2.0d0*DelTime &
& * wz_DiffVorDiv * wa_xya( xyz_Vor_b ) &
& )
!!$ call io_gt4_out_Put('xyz_Vor_Diffusion', &
!!$ & xya_wa( 2.0d0*DelTime &
!!$ & * wz_DiffVorDiv * wa_xya( xyz_Vor_b ) ) )
xyz_Div_a = &
& xya_wa( &
& wa_xya(xyz_Div_a) &
& + 2.0d0*DelTime &
& * wz_DiffVorDiv * wa_xya( xyz_Div_b ) &
& )
!!$ call io_gt4_out_Put('xyz_Div_Diffusion', &
!!$ & xya_wa( 2.0d0*DelTime &
!!$ & * wz_DiffVorDiv * wa_xya( xyz_Div_b ) ) )
xyz_Temp_a = &
& xya_wa( &
& wa_xya(xyz_Temp_a) &
& + 2.0d0*DelTime &
& * wz_DiffTherm * wa_xya( xyz_Temp_b ) &
& )
!!$ call io_gt4_out_Put('xyz_Temp_Diffusion', &
!!$ & xya_wa( 2.0d0*DelTime &
!!$ & * wz_DiffTherm * wa_xya( xyz_Temp_b ) ) )
xyz_QVap_a = &
& xya_wa( &
& wa_xya(xyz_QVap_a) &
& + 2.0d0*DelTime &
& * wz_DiffTherm * wa_xya( xyz_QVap_b ) &
& )
!!$ call io_gt4_out_Put('xyz_QVap_Diffusion', &
!!$ & xya_wa( 2.0d0*DelTime &
!!$ & * wz_DiffTherm * wa_xya( xyz_QVap_b ) ) )
! スペクトルデータとして見て拡散が効いているか確認するための
! データ出力を記述したソース。本計算にはまったく必要でない。
!!wz_Vor_a = wa_xya(xyz_Vor_a )
!!wz_Div_a = wa_xya(xyz_Div_a )
!!wz_Temp_a = wa_xya(xyz_Temp_a)
!!wz_QVap_a = wa_xya(xyz_QVap_a)
!!
!!write(80, *) CurrentTime, &
!! & wz_Vor_a( l_nm(21,0), 1) , wz_Div_a( l_nm(21,0), 1) , &
!! & wz_Temp_a( l_nm(21,0), 1) , wz_QVap_a( l_nm(21,0), 1)
!!call flush(80)
call EndSub(subname)
end subroutine
((<dynamics_init>)) で allocate した変数を deallocate し、 演算した値も全て破棄する。
subroutine dynamics_end ()
!==== Dependency
use type_mod, only: STRING, REKIND, DBKIND, INTKIND
use dc_trace, only: BeginSub, EndSub, DbgMessage
implicit none
!-----------------------------------------------------------------
! 変数定義
!-----------------------------------------------------------------
!----- 作業用内部変数 -----
character(STRING), parameter:: subname = "dynamics_end"
continue
!-----------------------------------------------------------------
! Check Initialization
!-----------------------------------------------------------------
call BeginSub(subname)
if ( .not. dynamics_initialized) then
call EndSub( subname, 'dynamics_init was not called', &
& c1=trim(subname) )
return
else
dynamics_initialized = .false.
endif
!-----------------------------------------------------------------
! allocate した変数の nullify とか
!-----------------------------------------------------------------
deallocate( xy_Coli , &
& xy_UVFact , &
& z_DelSigma , &
& wz_DiffVorDiv , &
& wz_DiffTherm , &
& z_HydroAlpha , &
& z_HydroBeta , &
& z_TempInpolA , &
& z_TempInpolB , &
& z_TempInpolKappa , &
& z_TempAve , &
& z_TempAveHalf , &
& zz_DivMtx , &
& zz_TempMtx , &
& zz_TempMtxQ , &
& zz_TempMtxS , &
& zz_TempMtxR )
deallocate &
& ( xy_lnPs , & ! π＝ ln Ｐs (t)
& xyz_lnPsAdv , & ! π の移流 ｖ・∇π
& xyz_lnPsAdvSum, & ! π移流積下げ Σ[k〜K](ｖ・∇π)Δσ
& xyz_DivSum , & ! 発散の積下げ Σ[k〜K] ＤΔσ
&
& xy_lnPs_T , & ! πの時間変化 dπ/dt (Δt)
& w_lnPs_T , & ! πの時間変化 (スペクトル) (Δt)
& w_lnPs_a , & ! π (t+Δt)
&
& xyz_VSigmaHalf , & ! 鉛直σ速度(半整数レベル) (t)
& xyz_VSigmaHalfNonG , & ! 鉛直σ速度 非重力波成分 (半整数レベル) (t)
&
& xyz_TempEdd , & ! Ｔ' ：温度の擾乱 (整数レベル) (t)
& xyz_TempVir , & ! Ｔv ：仮温度 (virtual temperature)(t)
& xyz_TempVirEdd, &! Ｔv'：仮温度の擾乱 (t)
& xyz_TempEddHalf, &! Ｔ' ：温度の擾乱 (半整数レベル) (t)
&
& xyz_UA_T , & ! 東西運動量変化項ＵＡ
& xyz_VA_T , & ! 南北運動量変化項ＶＡ
&
& wz_Vor_T , & ! 渦度変化のスペクトルデータ (Δt)
& wz_Vor_a , & ! 渦度のスペクトルデータ (t+Δt)
&
& xyz_KE , & ! 運動エネルギー + 仮温度補正
& ! ｕ**2+ｖ**2 + ΣＷ(Ｔv-Ｔ)
&
& wz_Div_T , & ! 発散のスペクトルデータ (Δt)
& wz_Div_a , & ! 発散のスペクトルデータ (t+Δt)
& wz_PresTendTemp , & ! 温度による圧力傾度のスペクトルデータ
& wz_PresTendPs , & ! 地表圧力による圧力傾度のスペクトルデータ
&
& xyz_TempLocal_T , & ! 温度の局所時間変化項 Ｈ
& ! (水平発散 Ｔ'Ｄ + σ移流
& ! + π変化の効果) (全て非重力波項)
&
& wz_Temp_T , & ! 温度の変化スペクトルデータ (Δt)
& wz_Temp_a , & ! 温度のスペクトルデータ (t+Δt)
&
& xyz_QVapDivVSigmaAdv, & ! Ｒ＝ｑＤ＋σ移流
& wz_QVap_T , & ! 比湿の変化スペクトルデータ (Δt)
& wz_QVap_a & ! 比湿のスペクトルデータ (t+Δt)
& )
deallocate &
& ( wz_Psi_a , & ! スペクトル(流線関数)
& wz_Chi_a & ! スペクトル(ポテンシャル)
& )
call EndSub(subname)
end subroutine
以降のサブルーチンで用いる変数の allocate 、および 時間発展しない量の演算を行なう。 また、変数データ出力のための初期設定も行なう。
subroutine dynamics_init (
x_Lon, y_Lat, z_Sigma, r_Sigma)
use type_mod, only: STRING, REKIND, DBKIND, INTKIND
use grid_3d_mod, only: im, jm, km
use grid_wavenumber_mod, only: nm
use constants_mod, only: constants_init, R0, Omega, Cp, RAir, &
& TempAve, VisOrder, EFoldTime
use time_mod, only: DelTime
use spml_mod, only: spml_init, xy_Lat, rn
use io_gt4_out_mod,only: io_gt4_out_init, io_gt4_out_SetVars
use dc_trace, only: DbgMessage, BeginSub, EndSub, DataDump
use dc_string, only: toChar
implicit none
real(DBKIND), intent(in) :: &
& x_Lon(:) , & ! intent(in): 経度座標
& y_Lat(:) , & ! intent(in): 緯度座標
& z_Sigma(:) , & ! intent(in): σレベル(整数)座標
& r_Sigma(:) ! intent(in): σレベル(半整数)座標
!----- 拡散係数計算用変数 -----
real(DBKIND) :: VisCoef ! 超粘性係数
real(DBKIND), pointer :: wz_rn(:,:) =>null() ! ラプラシアンの係数 -n*(n+1)
!----- 作業用内部変数 -----
! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
integer(INTKIND) :: i, j, k, kk, l
character(STRING), parameter:: subname = "dynamics_init"
continue
!----------------------------------------------------------------
! Check Initialization
!----------------------------------------------------------------
call BeginSub(subname)
if (dynamics_initialized) then
call EndSub( subname, '%c is already called.', c1=trim(subname) )
return
else
dynamics_initialized = .true.
endif
!----------------------------------------------------------------
! Version identifier
!----------------------------------------------------------------
call DbgMessage('%c :: %c', c1=trim(version), c2=trim(tagname))
!-------------------------------------------------------------------
! 物理定数の取得
!-------------------------------------------------------------------
call constants_init
!-------------------------------------------------------------------
! SPMODEL の初期化
!-------------------------------------------------------------------
call spml_init
!-------------------------------------------------------------------
! io_gt4_out の初期化
!-------------------------------------------------------------------
call io_gt4_out_init
!-------------------------------------------------------------------
! コリオリ力の設定
!-------------------------------------------------------------------
allocate( xy_Coli(im, jm) )
xy_Coli(:,:) = 2.0 * Omega * sin(xy_Lat(:,:))
!!$ call DataDump('xy_Coli', xy_Coli(:,:), strlen=80)
!!$ call DataDump('xy_Lat', xy_Lat(:,:), strlen=80)
!-------------------------------------------------------------------
! u→U, v→V のファクターの計算
!-------------------------------------------------------------------
allocate( xy_UVFact(im, jm) )
xy_UVFact(:,:) = cos( xy_Lat(:,:) )
!!$ call DataDump('xy_UVFact', xy_UVFact(:,:), strlen=80)
!-------------------------------------------------------------------
! Δσ (整数、半整数) の計算
!-------------------------------------------------------------------
allocate( z_DelSigma(km) )
do k = 1, km
z_DelSigma(k) = r_Sigma(k) - r_Sigma(k+1)
enddo
!!$ call DataDump('z_Delsigma', z_DelSigma(:), strlen=80)
!-------------------------------------------------------------------
! 運動量 (渦度発散) 拡散係数 wz_DiffVorDiv,
! 熱拡散係数 wz_DiffTherm の計算
!-------------------------------------------------------------------
allocate( wz_DiffVorDiv((nm+1)*(nm+1), km) )
allocate( wz_DiffTherm ((nm+1)*(nm+1), km) )
allocate( wz_rn ((nm+1)*(nm+1), km) )
do k = 1, km
wz_rn(:,k) = rn(:,1)
enddo
! 粘性係数の計算 (最大波数の e-folding time が EFoldTime となるように)
VisCoef = ( real( ( nm*(nm+1) ), DBKIND ) &
& / R0**2 )**(-VisOrder/2) &
& / EFoldTime
call DbgMessage('VisCoef=<%f>', d=(/VisCoef/))
wz_DiffTherm = &
& - VisCoef * ( ( - wz_rn / R0**2 )**(VisOrder/2) )
wz_DiffVorDiv = &
& wz_DiffTherm &
& - VisCoef &
& * ( - (2.0d0/R0**2)**(VisOrder/2) )
!!$ call DataDump('wz_DiffTherm' , wz_DiffTherm, strlen=80)
!!$ call DataDump('wz_DiffVorDiv', wz_DiffVorDiv, strlen=80)
deallocate(wz_rn)
!-------------------------------------------------------------------
! 静水圧の式の係数α、βの計算
!-------------------------------------------------------------------
allocate(z_HydroAlpha(km))
allocate(z_HydroBeta(km))
Kappa = RAir / Cp ! κ=Ｒ／Ｃｐ (気体定数/定圧比熱)
do k = 1, km
z_HydroAlpha(k) = &
& ( r_Sigma(k)/z_Sigma(k) )**Kappa - 1.0
z_HydroBeta(k) = &
& 1.0 - ( r_Sigma(k+1)/z_Sigma(k) )**Kappa
enddo
!!$ call DbgMessage('Kappa=<%f>', d=(/Kappa/) )
!!$ call DataDump('z_HydroAlpha', z_HydroAlpha(:), strlen=80)
!!$ call DataDump('z_HydroBeta', z_HydroBeta(:), strlen=80)
!-------------------------------------------------------------------
! 温度鉛直補間の係数ａ、ｂの計算
!-------------------------------------------------------------------
allocate(z_TempInpolA(km))
allocate(z_TempInpolB(km))
allocate(z_TempInpolKappa(km))
do k = 1, km
z_TempInpolA(k) = z_HydroAlpha(k) &
& / ( 1.0 - (z_Sigma(k) / z_Sigma(k-1))**Kappa )
z_TempInpolB(k) = z_HydroBeta(k) &
& / ( (z_Sigma(k) / z_Sigma(k+1))**Kappa - 1.0 )
z_TempInpolKappa(k) = &
& ( &
& r_Sigma(k) * z_HydroAlpha(k) &
& + r_Sigma(k+1) * z_HydroBeta(k) &
& ) / z_DelSigma(k)
enddo
z_TempInpolA(1) = 0.0
z_TempInpolB(km) = 0.0
!!$ call DataDump('z_TempInpolA', z_TempInpolA(:), strlen=80)
!!$ call DataDump('z_TempInpolB', z_TempInpolB(:), strlen=80)
!!$ call DataDump('z_TempInpolKappa', z_TempInpolKappa(:), strlen=80)
!-------------------------------------------------------------------
! 平均温度 (整数レベル、半整数レベル) の計算
!-------------------------------------------------------------------
allocate( z_TempAve(km) )
allocate( z_TempAveHalf(km+1) )
z_TempAve(:) = TempAve
!----- 下端の平均温度設定 -----
! ※ 正しいかは微妙だが一応変な値が入らないように
z_TempAveHalf(:) = 0.0d0
do k = 2, km
z_TempAveHalf(k) = &
& z_TempInpolA(k) * z_TempAve(k) &
& + z_TempInpolB(k-1) * z_TempAve(k-1)
enddo
!!$ call DataDump('z_TempAve', z_TempAve(:), strlen=80)
!!$ call DataDump('z_TempAveHalf', z_TempAveHalf(:), strlen=80)
!-------------------------------------------------------------------
! semi-implicit 用行列の計算
!-------------------------------------------------------------------
allocate( zz_DivMtx(km,km) )
allocate( zz_TempMtx(km,km) )
allocate( zz_TempMtxQ(km,km) )
allocate( zz_TempMtxS(km,km) )
allocate( zz_TempMtxR(km,km) )
! Ｗ:発散の式での温度補正 (線形重力波項の効果)
! 初期化
zz_DivMtx(:,:) = 0.0
do k = 1, km
do kk = 1, k
zz_DivMtx(k, kk) = Cp * z_HydroAlpha(kk)
enddo
do kk = 1, k-1
zz_DivMtx(k, kk) = zz_DivMtx(k, kk) + Cp * z_HydroBeta(kk)
enddo
enddo
! Ｓ: dσ/dｔ (線形重力波項の効果)
! 初期化
zz_TempMtxS(:,:) = 0.0
do k = 1, km
do kk = 1, km
zz_TempMtxS(k,kk) = r_Sigma(k) * z_DelSigma(kk)
enddo
do kk = k, km
zz_TempMtxS(k,kk) = zz_TempMtxS(k,kk) - z_DelSigma(kk)
enddo
enddo
!!$ call DataDump('zz_TempMtxS', zz_TempMtxS(:,:), strlen=80)
! Ｑ: dＴ/dσ (線形重力波項の効果)
! 初期化
zz_TempMtxQ(:,:) = 0.0
do k = 1, km
zz_TempMtxQ(k,k) = ( z_TempAveHalf(k) - z_TempAve(k) ) / z_DelSigma(k)
enddo
do k = 1, km - 1
zz_TempMtxQ(k,k+1) = ( z_TempAve(k) - z_TempAveHalf(k+1) ) / z_DelSigma(k)
enddo
!!$ call DataDump('zz_TempMtxQ', zz_TempMtxQ(:,:), strlen=80)
! Ｒ: ＲＤ = κＴ（∂π/∂ｔ＋ｄσ/ｄｔ／σ）
! 気圧変化の効果の係数 (線形重力波項の効果)
! 初期化
zz_TempMtxR(:,:) = 0.0
do k = 1, km
do kk = k, km
zz_TempMtxR(k,kk) = &
& - z_HydroAlpha(k) / z_DelSigma(k) &
& * z_DelSigma(kk) * z_TempAve(k)
enddo
do kk = k+1, km
zz_TempMtxR(k,kk) = zz_TempMtxR(k,kk) &
& - z_HydroBeta(k) / z_DelSigma(k) &
& * z_DelSigma(kk) * z_TempAve(k)
enddo
enddo
!!$ call DataDump('zz_TempMtxR', zz_TempMtxR(:,:), strlen=80)
! ｈ: ＱＳ (σ移流) - Ｒ
! 初期化
zz_TempMtx(:,:) = 0.0
zz_TempMtx(:,:) = &
& matmul(zz_TempMtxQ(:,:), zz_TempMtxS(:,:)) - zz_TempMtxR(:,:)
!!$ call DataDump('zz_TempMtx', zz_TempMtx(:,:), strlen=80)
!-------------------------------------------------------------------
! データ出力用の初期設定
!-------------------------------------------------------------------
call io_gt4_out_SetVars('xyz_DivSum')
call io_gt4_out_SetVars('xyz_lnPsAdv')
call io_gt4_out_SetVars('xyz_lnPsAdvSum')
call io_gt4_out_SetVars('xyz_UA_T')
call io_gt4_out_SetVars('xyz_VA_T')
call io_gt4_out_SetVars('xyz_KE')
call io_gt4_out_SetVars('xyz_PresTendTemp')
call io_gt4_out_SetVars('xyz_PresTendPs')
call io_gt4_out_SetVars('xyz_VSigmaHalf')
call io_gt4_out_SetVars('xyz_VSigmaHalfNonG')
call io_gt4_out_SetVars('xyz_TempLocal_T')
call io_gt4_out_SetVars('xyz_TempAdv')
call io_gt4_out_SetVars('xyz_TempEddHalf')
call io_gt4_out_SetVars('xyz_Temp_T')
call io_gt4_out_SetVars('xyz_TempLocal_T_lnPsAdvSum')
call io_gt4_out_SetVars('TotalMass')
call io_gt4_out_SetVars('xyz_Vor_Diffusion')
call io_gt4_out_SetVars('xyz_Div_Diffusion')
call io_gt4_out_SetVars('xyz_Temp_Diffusion')
call io_gt4_out_SetVars('xyz_QVap_Diffusion')
!-------------------------------------------------------------------
! dynamics_leapfrog で計算する変数の allocate
!-------------------------------------------------------------------
allocate &
& ( xy_lnPs (im, jm) , & ! π＝ ln Ｐs (t)
& xyz_lnPsAdv(im, jm, km) , & ! π の移流 ｖ・∇π
& xyz_lnPsAdvSum(im, jm, km), & ! π移流積下げ Σ[k〜K](ｖ・∇π)Δσ
& xyz_DivSum(im, jm, km) , & ! 発散の積下げ Σ[k〜K] ＤΔσ
&
& xy_lnPs_T( im, jm ) , & ! πの時間変化 dπ/dt (Δt)
& w_lnPs_T( (nm+1)*(nm+1) ) , & ! πの時間変化 (スペクトル) (Δt)
& w_lnPs_a( (nm+1)*(nm+1) ) , & ! π (t+Δt)
&
& xyz_VSigmaHalf(im, jm, km+1) , & ! 鉛直σ速度(半整数レベル) (t)
& xyz_VSigmaHalfNonG(im, jm, km+1) , & ! 鉛直σ速度 非重力波成分 (半整数レベル) (t)
&
& xyz_TempEdd(im, jm, km) , & ! Ｔ' ：温度の擾乱 (整数レベル) (t)
& xyz_TempVir(im, jm, km) , & ! Ｔv ：仮温度 (virtual temperature)(t)
& xyz_TempVirEdd(im, jm, km ), &! Ｔv'：仮温度の擾乱 (t)
& xyz_TempEddHalf(im, jm, km+1), &! Ｔ' ：温度の擾乱 (半整数レベル) (t)
&
& xyz_UA_T(im, jm, km) , & ! 東西運動量変化項ＵＡ
& xyz_VA_T(im, jm, km) , & ! 南北運動量変化項ＶＡ
&
& wz_Vor_T( (nm+1)*(nm+1), km ) , & ! 渦度変化のスペクトルデータ (Δt)
& wz_Vor_a( (nm+1)*(nm+1), km ) , & ! 渦度のスペクトルデータ (t+Δt)
&
& xyz_KE(im, jm, km) , & ! 運動エネルギー + 仮温度補正
& ! ｕ**2+ｖ**2 + ΣＷ(Ｔv-Ｔ)
&
& wz_Div_T( (nm+1)*(nm+1), km) , & ! 発散のスペクトルデータ (Δt)
& wz_Div_a( (nm+1)*(nm+1), km) , & ! 発散のスペクトルデータ (t+Δt)
& wz_PresTendTemp( (nm+1)*(nm+1), km) , & ! 温度による圧力傾度のスペクトルデータ
& wz_PresTendPs( (nm+1)*(nm+1), km) , & ! 地表圧力による圧力傾度のスペクトルデータ
&
& xyz_TempLocal_T(im, jm, km) , & ! 温度の局所時間変化項 Ｈ
& ! (水平発散 Ｔ'Ｄ + σ移流
& ! + π変化の効果) (全て非重力波項)
&
& wz_Temp_T( (nm+1)*(nm+1), km ) , & ! 温度の変化スペクトルデータ (Δt)
& wz_Temp_a( (nm+1)*(nm+1), km ) , & ! 温度のスペクトルデータ (t+Δt)
&
& xyz_QVapDivVSigmaAdv( im, jm, km), & ! Ｒ＝ｑＤ＋σ移流
& wz_QVap_T( (nm+1)*(nm+1), km) , & ! 比湿の変化スペクトルデータ (Δt)
& wz_QVap_a( (nm+1)*(nm+1), km) & ! 比湿のスペクトルデータ (t+Δt)
& )
!-------------------------------------------------------------------
! dynamics_diagnostic で計算する変数の allocate
!-------------------------------------------------------------------
allocate &
& ( wz_Psi_a( (nm+1)*(nm+1), km) , & ! スペクトル(流線関数)
& wz_Chi_a( (nm+1)*(nm+1), km) ) ! スペクトル(ポテンシャル)
call EndSub(subname)
end subroutine
時間発展する量の演算を行なう。 演算するデータの出力も行なう。
subroutine dynamics_leapfrog (
x_Lon , y_Lat , z_Sigma , r_Sigma ,
xyz_VelLon_b, xyz_VelLat_b, xyz_Vor_b, xyz_Div_b ,
xyz_Temp_b , xyz_QVap_b , xy_Ps_b ,
xyz_VelLon_n, xyz_VelLat_n, xyz_Vor_n, xyz_Div_n ,
xyz_Temp_n , xyz_QVap_n , xy_Ps_n ,
xyz_VelLon_a, xyz_VelLat_a, xyz_Vor_a, xyz_Div_a ,
xyz_Temp_a , xyz_QVap_a , xy_Ps_a )
!==== Dependency
use type_mod, only: STRING, REKIND, DBKIND, INTKIND
use grid_3d_mod, only: im, jm, km
use grid_wavenumber_mod, only: nm
use constants_mod, only: R0, Cp, EpsVT
use time_mod, only: DelTime, CurrentTime
use spml_mod, only: w_xy, xy_w , xy_GradLon_w, xy_GradLat_w, &
& w_Div_xy_xy, w_LaplaInv_w, &
& wa_xya, xya_wa, wa_Div_xya_xya, wa_Lapla_wa
use io_gt4_out_mod,only: io_gt4_out_Put
use dc_trace, only: DbgMessage, BeginSub, EndSub, DataDump
use dc_string, only: toChar
implicit none
!==== Input
!
real(DBKIND), intent(in) :: &
& x_Lon(:) , & ! intent(in): 経度座標
& y_Lat(:) , & ! intent(in): 緯度座標
& z_Sigma(:) , & ! intent(in): σレベル(整数)座標
& r_Sigma(:) , & ! intent(in): σレベル(半整数)座標
&
& xyz_VelLon_b(:,:,:) , & ! intent(in): 速度経度成分 (t-Δt)
& xyz_VelLat_b(:,:,:) , & ! intent(in): 速度緯度成分 (t-Δt)
& xyz_Vor_b(:,:,:) , & ! intent(in): 渦度 (t-Δt)
& xyz_Div_b(:,:,:) , & ! intent(in): 発散 (t-Δt)
& xyz_Temp_b(:,:,:) , & ! intent(in): 温度 (t-Δt)
& xyz_QVap_b(:,:,:) , & ! intent(in): 比湿 (t-Δt)
& xy_Ps_b(:,:) , & ! intent(in): 地表面気圧 (t-Δt)
&
& xyz_VelLon_n(:,:,:) , & ! intent(in): 速度経度成分 (t)
& xyz_VelLat_n(:,:,:) , & ! intent(in): 速度緯度成分 (t)
& xyz_Vor_n(:,:,:) , & ! intent(in): 渦度 (t)
& xyz_Div_n(:,:,:) , & ! intent(in): 発散 (t)
& xyz_Temp_n(:,:,:) , & ! intent(in): 温度 (t)
& xyz_QVap_n(:,:,:) , & ! intent(in): 比湿 (t)
& xy_Ps_n(:,:) ! intent(in): 地表面気圧 (t)
!==== Output
!
real(DBKIND), intent(out) :: &
& xyz_VelLon_a(:,:,:) , & ! intent(out): 速度経度成分 (t+Δt)
& xyz_VelLat_a(:,:,:) , & ! intent(out): 速度緯度成分 (t+Δt)
& xyz_Vor_a(:,:,:) , & ! intent(out): 渦度 (t+Δt)
& xyz_Div_a(:,:,:) , & ! intent(out): 発散 (t+Δt)
& xyz_Temp_a(:,:,:) , & ! intent(out): 温度 (t+Δt)
& xyz_QVap_a(:,:,:) , & ! intent(out): 比湿 (t+Δt)
& xy_Ps_a(:,:) ! intent(out): 地表面気圧 (t+Δt)
!----- 作業用内部変数 -----
! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
integer(INTKIND) :: i, j, k, kk, l
character(STRING), parameter:: subname = "dynamics_leapfrog"
!----------------------------------------------------------------
! Check Initialization
!----------------------------------------------------------------
call BeginSub(subname)
if (.not. dynamics_initialized) then
call EndSub( subname, 'Call dynamics_init before call %c', &
& c1=trim(subname) )
return
endif
!-------------------------------------------------------------------
! 地表気圧変化・鉛直σ速度の計算 (連続の式、熱の式で利用される)
!-------------------------------------------------------------------
xy_lnPs(:,:) = log( xy_Ps_n(:,:) )
! 地表面気圧 π の移流 = ｖ・∇π
do k = 1, km
xyz_lnPsAdv(:,:,k) = &
& ( &
& xyz_VelLon_n(:,:,k) &
& * xy_GradLon_w( w_xy(xy_lnPs(:,:)) ) &
&
& + xyz_VelLat_n(:,:,k) &
& * xy_GradLat_w( w_xy(xy_lnPs(:,:)) ) &
&
& ) / R0
enddo
call io_gt4_out_Put('xyz_lnPsAdv', xyz_lnPsAdv(:,:,:))
!!$ call DataDump('xyz_lnPsAdv', xyz_lnPsAdv(:,:,:), strlen=80)
! π移流の積み下げ Σ[k〜K] (ｖ・∇π) Δσ
xyz_lnPsAdvSum(:,:,km) = xyz_lnPsAdv(:,:, km ) * z_DelSigma( km )
do k = km-1 , 1, -1
xyz_lnPsAdvSum(:,:,k) = &
& xyz_lnPsAdvSum(:,:,k+1) &
& + xyz_lnPsAdv(:,:,k) * z_DelSigma(k)
enddo
call io_gt4_out_Put('xyz_lnPsAdvSum', xyz_lnPsAdvSum(:,:,:))
!!$ call DataDump('xyz_lnPsAdvSum', xyz_lnPsAdvSum(:,:,:), strlen=80)
! 発散の積下げ Σ[k〜K] ＤΔσ
xyz_DivSum(:,:,km) = &
& xyz_Div_n(:,:, km ) * z_DelSigma( km )
do k = km -1 , 1, -1
xyz_DivSum(:,:,k) = &
& xyz_DivSum(:,:,k+1) &
& + xyz_Div_n(:,:,k) * z_DelSigma(k)
enddo
call io_gt4_out_Put('xyz_DivSum', xyz_DivSum(:,:,:))
!!$ call DataDump('xyz_DivSum', xyz_DivSum(:,:,:), strlen=80)
!-------------------------------------------------------------------
! Ps の計算。連続の式を解く
!-------------------------------------------------------------------
!----- π = lnＰs の時間変化 dπ/dt の計算 -----
!
! π移流の積み下げ Σ[1〜K] (ｖ・∇π) Δσ を格子点上で解く。
!
xy_lnPs_T(:,:) = &
& - xyz_lnPsAdvSum(:,:,1) ! - Σ[1〜K](ｖ・∇π)Δσ
!!$ call DataDump('xy_lnPs_T', xy_lnPs_T(:,:), strlen=80)
!----- 発散の積み下げ -Σ[1〜K] ＤΔσ をスペクトル空間で加える。-----
w_lnPs_T(:) = w_xy( xy_lnPs_T(:,:) ) - w_xy( xyz_DivSum(:,:,1) )
!----- A = B + 2Δt * T-----
w_lnPs_a(:) = &
& w_xy( log( xy_Ps_b(:,:) ) ) &
& + 2.0*DelTime * ( w_lnPs_T(:) )
!!$ call DataDump('xy_Ps_b', xy_Ps_b(:,:), strlen=80)
!!$ call DataDump('w_lnPs_b', w_xy( log( xy_Ps_b(:,:) ) ), strlen=80)
!!$ call DataDump('w_lnPs_a', w_lnPs_a(:), strlen=80)
!----- xy_Ps_a に代入-----
xy_Ps_a(:,:) = exp( xy_w( w_lnPs_a(:) ) )
!!$ call DataDump('xy_Ps_a', xy_Ps_a(:,:), strlen=80)
!-------------------------------------------------------------------
! 鉛直σ速度 (整数、半整数レベル) の上昇流の計算。
!-------------------------------------------------------------------
! σ速度 (3/2 〜 K-1/2)
do k = 2, km+1 - 1
xyz_VSigmaHalf(:,:,k) = &
&
! - σ[k-1/2] × ( Σ[1〜K](Ｄ + ｖ・∇π)Δσ )
& r_Sigma(k) &
& * ( xyz_lnPsAdvSum(:,:,1) + xyz_DivSum(:,:,1) ) &
&
! - Σ[k〜K] (Ｄ + ｖ・∇π) Δσ
& - ( xyz_lnPsAdvSum(:,:,k) + xyz_DivSum(:,:,k) )
xyz_VSigmaHalfNonG(:,:,k) = &
&
! - σ[k-1/2] × ( Σ[1〜K](ｖ・∇π)Δσ )
& r_Sigma(k) * xyz_lnPsAdvSum(:,:,1) &
&
! - Σ[k〜K] (ｖ・∇π) Δσ
& - xyz_lnPsAdvSum(:,:,k)
enddo
! σ速度 (1/2, K+1/2)
xyz_VSigmaHalf(:,:,1) = 0.0
xyz_VSigmaHalf(:,:,km+1) = 0.0
xyz_VSigmaHalfNonG(:,:,1) = 0.0
xyz_VSigmaHalfNonG(:,:,km+1) = 0.0
call io_gt4_out_Put('xyz_VSigmaHalf', xyz_VSigmaHalf(:,:,:))
call io_gt4_out_Put('xyz_VSigmaHalfNonG', xyz_VSigmaHalfNonG(:,:,:))
!!$ call DataDump('xyz_VSigmaHalf', xyz_VSigmaHalf(:,:,:), strlen=80)
!-------------------------------------------------------------------
! 温度・仮温度の基本場からのずれ (以降の ＵＡ、ＵＶ、Ｈなどの全てで利用)
!-------------------------------------------------------------------
do k = 1, km
xyz_TempVir(:,:,k) = xyz_Temp_n(:,:,k) &
& * ( 1.0 + (EpsVT * xyz_QVap_n(:,:,k)) )
xyz_TempEdd(:,:,k) = xyz_Temp_n(:,:,k) - z_TempAve(k)
xyz_TempVirEdd(:,:,k) = xyz_TempVir(:,:,k) - z_TempAve(k)
enddo
!!$ call DataDump('xyz_Temp', xyz_Temp(:,:,:), strlen=80)
!!$ call DataDump('xyz_TempEdd', xyz_TempEdd(:,:,:), strlen=80)
!!$ call DataDump('xyz_TempVir', xyz_TempVir(:,:,:), strlen=80)
!!$ call DataDump('xyz_TempVirEdd', xyz_TempVirEdd(:,:,:), strlen=75)
!-------------------------------------------------------------------
! 半整数レベルの温度の擾乱 (以降の ＵＡ、ＵＶ、Ｈなどの全てで利用)
!-------------------------------------------------------------------
do k = 2, km+1 - 1
xyz_TempEddHalf(:,:,k) = &
& z_TempInpolA(k) * xyz_Temp_n(:,:,k) &
& + z_TempInpolB(k-1) * xyz_Temp_n(:,:,k-1) &
& - z_TempAveHalf(k)
enddo
call io_gt4_out_Put('xyz_TempEddHalf', xyz_TempEddHalf(:,:,:))
!!$ call DataDump('xyz_TempEddHalf', xyz_TempEddHalf(:,:,:), strlen=80)
!-------------------------------------------------------------------
! ＵＡ＝Ｖζ＋σ移流、ＶＡ＝Ｕζ＋σ移流 (渦度と発散の式で利用)
! ※ AGCM5 のマニュアルの ＵＡ と異なり、UVFact = cosφは
! 掛かっていないので注意。
!-------------------------------------------------------------------
do k = 1, km
xyz_UA_T(:,:,k) = &
! (ζ + ｆ) ｖ
& ( xyz_Vor_n(:,:,k) + xy_Coli(:,:) ) * xyz_VelLat_n(:,:,k) &
&
! - (Ｃp κ (1/ａ) Ｔ'v (dπ/dλ) ) / cosφ
& - Cp * z_TempInpolKappa(k) &
& * xyz_TempVirEdd(:,:,k) &
& * xy_GradLon_w( w_xy( xy_lnPs(:,:) ) ) &
& / R0
xyz_VA_T(:,:,k) = &
! - (ζ + ｆ) ｕ
& - ( xyz_Vor_n(:,:,k) + xy_Coli(:,:) ) * xyz_VelLon_n(:,:,k) &
&
! - (Ｃp κ (1/ａ) Ｔ'v (dπ/dμ) ) / cosφ
& - Cp * z_TempInpolKappa(k) &
& * xyz_TempVirEdd(:,:,k) &
& * xy_GradLat_w( w_xy( xy_lnPs(:,:) ) ) &
& / R0
enddo
!!$ do k = 1, km
!!$ call DataDump('xyz_Div(pi)', &
!!$ & xy_w( &
!!$ & w_Div_xy_xy( &
!!$ & - Cp * z_TempInpolKappa(k) &
!!$ & * xyz_TempVirEdd(:,:,k) &
!!$ & * xy_GradLon_w( w_xy( xy_lnPs(:,:) ) ) &
!!$ & / R0 &
!!$ & , &
!!$ &
!!$ & - Cp * z_TempInpolKappa(k) &
!!$ & * xyz_TempVirEdd(:,:,k) &
!!$ & * xy_GradLat_w( w_xy( xy_lnPs(:,:) ) ) &
!!$ & / R0 &
!!$ &
!!$ & ) &
!!$ & ) / R0 &
!!$ &
!!$ & , strlen=80, multi=(/k/))
!!$ enddo
!!$ call DataDump('xyz_VelLon', xyz_VelLon(:,:,:), strlen=80)
!!$ call DataDump('xyz_VelLat', xyz_VelLat(:,:,:), strlen=80)
!!$ call DataDump('xyz_VA_T', xyz_VA_T(1,1,:), strlen=80, multi=(/-1/))
!!$ call DataDump('xyz_UA_T', xyz_UA_T(1,1,:), strlen=80, multi=(/-1/))
!!$ call DataDump('xyz_VA_T', xyz_VA_T(1,1,:), strlen=80, multi=(/-1/))
do k = 2, km
xyz_UA_T(:,:,k) = xyz_UA_T(:,:,k) &
! - (1/2Δσ[k]) * (dσ/dt)[k-1/2] * (ｕ[k-1] - ｕ[k])
& - 1.0 / (2.0 * z_DelSigma(k)) &
& * xyz_VSigmaHalf(:,:,k) &
& * ( xyz_VelLon_n(:,:,k-1) - xyz_VelLon_n(:,:,k) )
xyz_VA_T(:,:,k) = xyz_VA_T(:,:,k) &
! - (1/2Δσ[k]) * (dσ/dt)[k-1/2] * (ｖ[k-1] - ｖ[k])
& - 1.0 / (2.0 * z_DelSigma(k)) &
& * xyz_VSigmaHalf(:,:,k) &
& * ( xyz_VelLat_n(:,:,k-1) - xyz_VelLat_n(:,:,k) )
enddo
!!$ call DataDump('xyz_UA_T', xyz_UA_T(1,1,:), strlen=80, multi=(/-2/))
!!$ call DataDump('xyz_VA_T', xyz_VA_T(1,1,:), strlen=80, multi=(/-2/))
do k = 1, km - 1
xyz_UA_T(:,:,k) = xyz_UA_T(:,:,k) &
! - (1/2Δσ[k]) * (dσ/dt)[k+1/2] * (ｕ[k] - ｕ[k+1])
& - 1.0 / (2.0 * z_DelSigma(k)) &
& * xyz_VSigmaHalf(:,:,k+1) &
& * ( xyz_VelLon_n(:,:,k) - xyz_VelLon_n(:,:,k+1) )
xyz_VA_T(:,:,k) = xyz_VA_T(:,:,k) &
! - (1/2Δσ[k]) * (dσ/dt)[k+1/2] * (ｖ[k] - ｖ[k+1])
& - 1.0 / (2.0 * z_DelSigma(k)) &
& * xyz_VSigmaHalf(:,:,k+1) &
& * ( xyz_VelLat_n(:,:,k) - xyz_VelLat_n(:,:,k+1) )
enddo
!!$ call DataDump('xyz_UA_T', xyz_UA_T(1,1,:), strlen=80, multi=(/-3/))
!!$ call DataDump('xyz_VA_T', xyz_VA_T(1,1,:), strlen=80, multi=(/-3/))
!!$ call DataDump('xyz_UA_T', xyz_UA_T(:,:,:), strlen=80)
!!$ call DataDump('xyz_VA_T', xyz_VA_T(:,:,:), strlen=80)
call io_gt4_out_Put('xyz_UA_T', xyz_UA_T(:,:,:))
call io_gt4_out_Put('xyz_VA_T', xyz_VA_T(:,:,:))
!-------------------------------------------------------------------
! 渦度 Vor の計算。渦度方程式を解く
!-------------------------------------------------------------------
!!$ ! dζ/dt の計算
!!$ wz_VorTend_n = &
!!$ & wa_Div_xya_xya( xyz_VA_n, - xyz_UA_n ) / Rplanet
!!$ !& + w_DiffV(k) * wa_xya(xyz_Vor(:,:,k))
!!$ ! 拡散係数の導出が厄介なので後回し
! dζ/dt の計算
wz_Vor_T(:,:) = &
& wa_Div_xya_xya( xyz_VA_T(:,:,:), - xyz_UA_T(:,:,:) ) / R0
!& + w_DiffV(k) * wa_xya(xyz_Vor(:,:,k))
! 拡散係数の導出が厄介なので後回し
!----- A = B + 2Δt * T-----
wz_Vor_a(:,:) = &
& wa_xya( xyz_Vor_b(:,:,:) ) &
& + 2.0d0*DelTime * ( wz_Vor_T(:,:) )
!!! & + 2.0d0*DelTime * wa_NumVis_wa( wa_xya( xyz_Vor_b(:,:,:) ) )
!----- xyz_Vor_a に代入-----
xyz_Vor_a(:,:,:) = xya_wa( wz_Vor_a(:,:) )
!!$ call DataDump('wz_Vor_T', wz_Vor_T(:,:), strlen=70)
!!$ call DataDump('xyz_Vor_b', xyz_Vor_b(:,:,:), strlen=80)
!!$ call DataDump('xyz_Vor', xyz_Vor(:,:,:), strlen=80)
!!$ call DataDump('xyz_Vor_T', xya_wa(wz_Vor_T(:,:)), strlen=80)
!!$ call DataDump('xyz_Vor_a', xyz_Vor_a(:,:,:), strlen=80)
!-------------------------------------------------------------------
! Ｅ=ｕ**2+ｖ**2 + 仮温度補正 ( ΣＷ(Ｔv-Ｔ) ) (発散の式で利用)
!-------------------------------------------------------------------
! 仮温度補正 ( ΣＷ(Ｔv-Ｔ) ) の計算。格子点上で積み上げ。
! 初期化
xyz_KE(:,:,:) = 0.0
xyz_KE(:,:,1) = Cp * z_HydroAlpha(1) &
& * ( xyz_TempVir(:,:,1) - xyz_Temp_n(:,:,1) )
do k = 2, km
xyz_KE(:,:,k) = xyz_KE(:,:,k-1) &
& + Cp * z_HydroAlpha(k) &
& * ( xyz_TempVir(:,:,k) - xyz_Temp_n(:,:,k) ) &
& + Cp * z_HydroBeta(k-1) &
& * ( xyz_TempVir(:,:,k-1) - xyz_Temp_n(:,:,k-1) )
enddo
!!$ call DataDump('xyz_Temp', xyz_Temp(:,:,:), strlen=80)
!!$ call DataDump('xyz_TempVir', xyz_TempVir(:,:,:), strlen=80)
!!$ call DataDump('xyz_KE-w(Tv-T)', xyz_KE(:,:,:), strlen=80)
! 運動エネルギー ｕ**2+ｖ**2 の加算
xyz_KE(:,:,:) = xyz_KE(:,:,:) &
& + ( &
& xyz_VelLon_n(:,:,:)**2 &
& + xyz_VelLat_n(:,:,:)**2 &
& ) / 2.0
!!$ call DataDump('xyz_KE-u**2+v**2', xyz_KE(:,:,:), strlen=80)
call io_gt4_out_Put('xyz_KE', xyz_KE(:,:,:))
!-------------------------------------------------------------------
! 発散 Div の計算。発散方程式を解く
!-------------------------------------------------------------------
! ＷＴ = Ｃp α Ｔ[k] + Ｃp β Ｔ[k-1] をスペクトル空間で積み上げ
! 初期化
wz_PresTendTemp(:,:) = 0.0
wz_PresTendTemp(:,1) = Cp * z_HydroAlpha(1) * w_xy( xyz_Temp_n(:,:,1) )
do k = 2, km
wz_PresTendTemp(:,k) = wz_PresTendTemp(:,k-1) &
& + Cp * z_HydroAlpha(k) * w_xy( xyz_Temp_n(:,:,k) ) &
& + Cp * z_HydroBeta(k-1) * w_xy( xyz_Temp_n(:,:,k-1) )
enddo
!!$ call io_gt4_out_Put('xyz_PresTendTemp', xya_wa(wa_Lapla_wa(wz_PresTendTemp(:,:)))/R0**2)
!!$ call DataDump('wz_PresTendTemp', wz_PresTendTemp(:,:), strlen=70)
!!$ call DataDump('xyz_Lapla(wz_PresTendTemp)/R0**2', &
!!$ & xya_wa( wa_Lapla_wa( wz_PresTendTemp(:,:) ) )/R0**2, strlen=60)
! Ｇπ = (Ｃp κ Ｔ) π
do k = 1, km
wz_PresTendPs(:,k) = &
& Cp * z_TempInpolKappa(k) * z_TempAve(k) &
& * w_xy( log( xy_Ps_n(:,:) ) )
enddo
!!$ call io_gt4_out_Put('xyz_PresTendPs', xya_wa(wa_Lapla_wa(wz_PresTendPs(:,:)))/R0**2)
!!$ call DataDump('wz_PresTendPs', wz_PresTendPs(:,:), strlen=70)
!!$ call DataDump('xyz_Lapla(wz_PresTendPs)/R0**2', &
!!$ & xya_wa( wa_Lapla_wa( wz_PresTendPs(:,:) ) )/R0**2, strlen=60)
!!$ call DataDump('Div(wz_UA_T - wz_VA_T)', &
!!$ & wa_Div_xya_xya( xyz_UA_T(:,:,:), - xyz_VA_T(:,:,:) ), &
!!$ & strlen=70)
!!$ call DataDump('Div(wz_UA_T - wz_VA_T)/R0', &
!!$ & wa_Div_xya_xya( xyz_UA_T(:,:,:), - xyz_VA_T(:,:,:) ) /R0, &
!!$ & strlen=70)
!!$ call DataDump('xyz_Div(wz_UA_T - wz_VA_T)/R0', &
!!$ & xya_wa( wa_Div_xya_xya( xyz_UA_T(:,:,:), - xyz_VA_T(:,:,:) ) ) /R0, &
!!$ & strlen=70)
!!$ call DataDump('Lapla(wz_KE)', &
!!$ & wa_Lapla_wa( wa_xya(xyz_KE(:,:,:)) ), strlen=70)
! dＤ/dt の計算
wz_Div_T(:,:) = &
& wa_Div_xya_xya( xyz_UA_T(:,:,:), xyz_VA_T(:,:,:) ) / R0 &
& - wa_Lapla_wa( wa_xya(xyz_KE(:,:,:)) ) / R0**2 &
!& + w_DiffV(k) * wa_xya(xyz_Div(:,:,k))
! 拡散係数の導出が厄介なので後回し
&
! semi-implicitだとここには無い項達
! - ∇**2 (Φ + ＷＴ + Ｇπ)
& - wa_Lapla_wa( &
!& w_xy( xy_Phi(:,:) ) & ! 地表ジオポテンシャルΦは後回し
& wz_PresTendTemp(:,:) & ! ＷＴ
& + wz_PresTendPs(:,:) & ! Ｇπ = (Ｃp κ Ｔ)π
& ) / R0**2
!----- A = B + 2Δt * T-----
wz_Div_a(:,:) = &
& wa_xya( xyz_Div_b(:,:,:) ) &
& + 2.0d0*DelTime * ( wz_Div_T(:,:) )
!!! & + 2.0d0*DelTime * wa_NumVis_wa( wa_xya( xyz_Div_b(:,:,:) ) )
!----- xyz_Div_a に代入-----
xyz_Div_a(:,:,:) = xya_wa( wz_Div_a(:,:) )
!!$ call DataDump('xyz_Div_b', xyz_Div_b(:,:,:), strlen=80)
!!$ call DataDump('xyz_Div', xyz_Div(:,:,:), strlen=80)
!!$ call DataDump('xyz_Div_T', xya_wa(wz_Div_T(:,:)), strlen=80)
!!$ call DataDump('xyz_Div_a', xyz_Div_a(:,:,:), strlen=80)
!-------------------------------------------------------------------
! Ｈ＝水平発散 Ｔ'Ｄ + σ移流 + π変化の効果 (熱力学の式)
!-------------------------------------------------------------------
! 温度の局所時間変化項 Ｈ の計算
do k = 1, km
xyz_TempLocal_T(:,:,k) = &
&
! Ｔ' Ｄ
& xyz_TempEdd(:,:,k) * xyz_Div_n(:,:,k) &
&
! κＴv (ｖ・∇π) : πの移流の効果
& + z_TempInpolKappa(k) * xyz_TempVir(:,:,k) &
& * xyz_lnPsAdv(:,:,k) &
&
! - (α/Δσ) {Ｔv Σ[k〜K] (ｖ・∇π) + Ｔv' Σ[k〜K] (ＤΔσ)}
& - z_HydroAlpha(k) / z_DelSigma(k) &
& * ( &
& xyz_TempVir(:,:,k) * xyz_lnPsAdvSum(:,:,k) &
& + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k) &
& )
enddo
do k = 2, km
xyz_TempLocal_T(:,:,k) = xyz_TempLocal_T(:,:,k) &
&
! - (dσ/dt) * (１/Δσ) * (Ｔ'[k-1/2] − Ｔ'[k])
& - xyz_VSigmaHalf(:,:,k) / z_DelSigma(k) &
& * ( xyz_TempEddHalf(:,:,k) - xyz_TempEdd(:,:,k) ) &
&
! - (dσ/dt)^NG * (１/Δσ) * (Ｔ￣[k-1/2] − Ｔ￣[k])
& - xyz_VSigmaHalfNonG(:,:,k) / z_DelSigma(k) &
& * ( z_TempAveHalf(k) - z_TempAve(k) )
enddo
do k = 1, km - 1
xyz_TempLocal_T(:,:,k) = xyz_TempLocal_T(:,:,k) &
&
! - (dσ/dt) * (１/Δσ) * (Ｔ'[k] − Ｔ'[k+1/2])
& - xyz_VSigmaHalf(:,:,k+1) / z_DelSigma(k) &
& * ( xyz_TempEdd(:,:,k) - xyz_TempEddHalf(:,:,k+1) ) &
&
! - (dσ/dt)^NG * (１/Δσ) * (Ｔ￣[k] − Ｔ￣[k+1/2])
& - xyz_VSigmaHalfNonG(:,:,k+1) / z_DelSigma(k) &
& * ( z_TempAve(k) - z_TempAveHalf(k+1) ) &
&
! - (β/Δσ) { Ｔv Σ[k+1〜K] (ｖ・∇π)
! + Ｔv' Σ[k+1〜K] (ＤΔσ)}
& - z_HydroBeta(k) / z_DelSigma(k) &
& * ( &
& xyz_TempVir(:,:,k) * xyz_lnPsAdvSum(:,:,k+1) &
& + xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k+1) &
& )
enddo
!!$ call io_gt4_out_Put('xyz_TempLocal_T_lnPsAdvSum', &
!!$ & - z_HydroAlpha(2) / z_DelSigma(2) &
!!$ & * ( &
!!$ & xyz_TempVir(:,:,2) * xyz_lnPsAdvSum(:,:,2) &
!!$ & ) &
!!$ & - z_HydroBeta(2) / z_DelSigma(2) &
!!$ & * ( &
!!$ & xyz_TempVir(:,:,2) * xyz_lnPsAdvSum(:,:,2+1) &
!!$ & ) &
!!$ & )
!!$ call DataDump('xyz_DivSum*TempVirEdd*(Alpha+Beta)', &
!!$ & - z_HydroAlpha(k) / z_DelSigma(k) &
!!$ & * ( xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k) ) &
!!$ & - z_HydroBeta(k) / z_DelSigma(k) &
!!$ & * ( xyz_TempVirEdd(:,:,k) * xyz_DivSum(:,:,k+1) ) &
!!$ & , strlen=70)
!!$ call DataDump('xyz_TempLocal_T', xyz_TempLocal_T(:,:,:), strlen=80)
call io_gt4_out_Put('xyz_TempLocal_T', xyz_TempLocal_T(:,:,:))
!-------------------------------------------------------------------
! ＵＴ' 、ＶＴ' (熱力学の式)
!-------------------------------------------------------------------
!!$ DO 5100 K = 1, KMAX
!!$ DO 5100 IJ = 1, IDIM*JDIM
!!$ GTUT( IJ,K ) = GAU ( IJ,K ) * GATED ( IJ,K )
!!$ & * UVFACT ( IJ )
!!$ GTVT( IJ,K ) = GAV ( IJ,K ) * GATED ( IJ,K )
!!$ & * UVFACT ( IJ )
!!$ 5100 CONTINUE
!-------------------------------------------------------------------
! 温度 Ｔ の計算。熱力学の式を解く
!-------------------------------------------------------------------
! 非重力波項成分 (非線形項) の部分に関する温度変化を計算
do k = 1, km
wz_Temp_T(:,k) = &
&
! dＵＴ'/dλ + dＶＴ'/dμ
& - w_Div_xy_xy( &
& xyz_VelLon_n(:,:,k) * xyz_TempEdd(:,:,k), &
& xyz_VelLat_n(:,:,k) * xyz_TempEdd(:,:,k) &
& ) / R0 &
&
! 温度の局所時間変化項 Ｈ
& + w_xy( xyz_TempLocal_T(:,:,k) )
!&
! 非断熱加熱 Ｑ もまだ入ってないよ
!&
! 摩擦熱および熱拡散は後回し
enddo
!!$ call DataDump('xyz_Temp_T_NonGrav', xya_wa(wz_Temp_T(:,:)), strlen=70)
! 重力波項成分 (線形項) の部分に関する温度変化を計算
do k = 1, km
do kk = 1, km
wz_Temp_T(:,k) = wz_Temp_T(:,k) &
&
! ｈＤ (重力波項成分によるπ変化の効果)
& - zz_TempMtx(k,kk) * w_xy( xyz_Div_n(:,:,kk) )
!&
! 拡散項は後回し
enddo
enddo
!----- A = B + 2Δt * T-----
wz_Temp_a(:,:) = &
& wa_xya( xyz_Temp_b(:,:,:) ) &
& + 2.0d0*DelTime * ( wz_Temp_T(:,:) )
!!! & + 2.0d0*DelTime * wa_NumVisScaler_wa( wa_xya( xyz_Temp_b(:,:,:) ) )
!----- xyz_Temp_a に代入-----
xyz_Temp_a(:,:,:) = xya_wa( wz_Temp_a(:,:) )
call io_gt4_out_Put( 'xyz_Temp_T', xya_wa( wz_Temp_T(:,:) ) )
!!$ call DataDump('xyz_Temp_b', xyz_Temp_b(:,:,:), strlen=80)
!!$ call DataDump('xyz_Temp', xyz_Temp(:,:,:), strlen=80)
!!$ call DataDump('xyz_Temp_T', xya_wa(wz_Temp_T(:,:)), strlen=80)
!!$ call DataDump('xyz_Temp_a', xyz_Temp_a(:,:,:), strlen=80)
!-------------------------------------------------------------------
! 比湿 QVap の計算。水蒸気の式を解く。
!-------------------------------------------------------------------
! Ｒ＝ｑＤ＋σ移流 (非重力波項) を格子点上で計算
xyz_QVapDivVSigmaAdv(:,:,:) = xyz_QVap_n(:,:,:) * xyz_Div_n(:,:,:)
do k = 2, km
xyz_QVapDivVSigmaAdv(:,:,k) = xyz_QVapDivVSigmaAdv(:,:,k) &
& - xyz_VSigmaHalf(:,:,k) / ( 2.0 * z_DelSigma(k) ) &
& * ( xyz_QVap_n(:,:,k-1) - xyz_QVap_n(:,:,k) )
enddo
do k = 1, km - 1
xyz_QVapDivVSigmaAdv(:,:,k) = xyz_QVapDivVSigmaAdv(:,:,k) &
& - xyz_VSigmaHalf(:,:,k+1) / ( 2.0 * z_DelSigma(k) ) &
& * ( xyz_QVap_n(:,:,k) - xyz_QVap_n(:,:,k+1) )
enddo
! 水平移流 dＵｑ/dλ + dＶｑ/dμ と水平拡散項、水蒸気ソースを計算
do k = 1, km
wz_QVap_T(:,k) = &
&
! dＵｑ/dλ + dＶｑ/dμ
& - w_Div_xy_xy( &
& xyz_VelLon_n(:,:,k) * xyz_QVap_n(:,:,k), &
& xyz_VelLat_n(:,:,k) * xyz_QVap_n(:,:,k) &
& ) / R0 &
&
! Ｒ＝ｑＤ＋σ移流
& + w_xy( xyz_QVapDivVSigmaAdv(:,:,k) )
!&
! 水平拡散項、水蒸気ソースは後回し
enddo
!----- A = B + 2Δt * T-----
wz_QVap_a(:,:) = &
& wa_xya( xyz_QVap_b(:,:,:) ) &
& + 2.0d0*DelTime * ( wz_QVap_T(:,:) )
!!! & + 2.0d0*DelTime * wa_NumVisScaler_wa( wa_xya( xyz_QVap_b(:,:,:) ) )
!----- xyz_QVap_a に代入-----
xyz_QVap_a(:,:,:) = xya_wa( wz_QVap_a(:,:) )
!!$ call DataDump('xyz_QVap_b', xyz_QVap_b(:,:,:), strlen=80)
!!$ call DataDump('xyz_QVap', xyz_QVap(:,:,:), strlen=80)
!!$ call DataDump('xyz_QVap_T', xya_wa(wz_QVap_T(:,:)), strlen=80)
!!$ call DataDump('xyz_QVap_a', xyz_QVap_a(:,:,:), strlen=80)
call EndSub(subname)
end subroutine