subroutine physics_cumulus_adjust( xyz_Temp , xyz_Qvap , xy_CumulusRain , xyz_DCumulusTempDt, xyz_DCumulusQvapDt, xyz_Press , xyr_Press , DelTimePhy ) ! (in) 2Δt
!==== Dependency
use type_mod, only: REKIND, DBKIND, INTKIND, TOKEN, STRING
use grid_3d_mod, only: im, jm, km
use constants_mod, only: Cp , EL , EpsV , ES0 , RVap , RAir , Grav ! 重力加速度
use physics_qvapsat_tetens, only: QvapSat
use dc_trace, only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump
implicit none
!==== Input
!
real(DBKIND), intent(in) :: xyz_Press(im,jm,km) , xyr_Press(im,jm,km+1) , DelTimePhy ! (in) 2Δt
!==== Output
!
real(DBKIND), intent(out) :: xy_CumulusRain(im,jm) , xyz_DCumulusTempDt(im,jm,km), xyz_DCumulusQvapDt(im,jm,km) ! (out) 比湿変化率
!==== In/Out
!
real(DBKIND), intent(inout) :: xyz_Temp(im,jm,km) , xyz_Qvap(im,jm,km) ! (inout) 比湿
!----- 作業用内部変数 -----
character(STRING), parameter:: subname = "physics_cumulus_adjust"
real(DBKIND) :: xyz_Qvap_b(im,jm,km) , xyz_Temp_b(im,jm,km) , xy_Adjust(im,jm) , xy_Adjust_b(im,jm) , xyz_DPressDz(im,jm,km) , xyz_QvapSat(im,jm,km) , xyz_DDPressDDPress(im,jm,km) , xyz_DPFact(im,jm,km) , TempSat , DelTempSat , DelQvap , DelTempUpper , DelTempUnder , DQvapSatDTempUpper , DQvapSatDTempUnder , DHDTempUpper , DHDTempUnder , Adjust ! 今回全領域において一度でも調節されたか否か
real(DBKIND), parameter :: CrtlRH = 0.990d0 ! 臨界相対湿度
integer(INTKIND), parameter :: IterationMax = 10 ! イテレーション回数
real(DBKIND) :: TempSatMax(IterationMax) ! 不安定の許容誤差
data TempSatMax / 0.01 , 0.02 , 0.02 , 0.05 , 0.05 , 0.10 , 0.10 , 0.20 , 0.20 , 0.40 /
! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
integer(INTKIND) :: i, j, k
integer(INTKIND) :: Iteration
continue
!----------------------------------------------------------------
! 開始処理
!----------------------------------------------------------------
call BeginSub(subname)
!----------------------------------------------------------------
! 湿潤対流調節
!----------------------------------------------------------------
!----- 調節前 Qvap の保存 -----
xyz_Qvap_b = xyz_Qvap
xyz_Temp_b = xyz_Temp
!----- 変数初期化 -----
xy_CumulusRain = 0.0d0
xyz_DCumulusTempDt = 0.0d0
xyz_DCumulusQvapDt = 0.0d0
!----- ファクター計算 -----
do k = 1, km
xyz_DPressDz(:,:,k) = xyr_Press(:,:,k) - xyr_Press(:,:,k+1)
end do
! クラウジウスクラペイロンの式より飽和比湿計算
! xyz_QvapSat = EpsV * ES0 &
! & * EXP( EL / RVap * ( 1./273. - 1./xyz_Temp ) ) / xyz_Press
do i = 1, im
do j = 1, jm
do k = 1, km
xyz_QvapSat(i,j,k) = QvapSat(xyz_Temp(i,j,k), xyz_Press(i,j,k))
enddo
enddo
enddo
do k = 2, km
xyz_DDPressDDPress(:,:,k) = xyz_DPressDz(:,:,k) / xyz_DPressDz(:,:,k-1)
xyz_DPFact(:,:,k) = RAir / Cp * ( xyz_Press(:,:,k-1) - xyz_Press(:,:,k) ) / xyr_Press(:,:,k) / 2.
end do
!----- 調節 ------
xy_Adjust_b = 1.0d0
!----- 3.1 イテレーション ------
do Iteration = 1, IterationMax
xy_Adjust = 0.0d0
do k = 2, km
do i = 1, im
do j = 1, jm
if ( xy_Adjust_b(i,j) .GT. 0.9 ) then
TempSat = xyz_Temp(i,j,k-1) - xyz_Temp(i,j,k) + ( xyz_QvapSat(i,j,k-1) - xyz_QvapSat(i,j,k) ) * EL/Cp - xyz_DPFact(i,j,k) * ( xyz_Temp(i,j,k-1) + xyz_Temp(i,j,k) )
! 不安定であるならば
if ( TempSat .GT. TempSatMax(Iteration) ) then
! かつ, 飽和しているならば
if ( ( xyz_Qvap(i,j,k)/xyz_QvapSat(i,j,k) .GE. CrtlRH ) .AND. ( xyz_Qvap(i,j,k-1)/xyz_QvapSat(i,j,k-1) .GE. CrtlRH ) ) then
DelQvap = xyz_DPressDz(i,j,k-1) * (xyz_Qvap(i,j,k-1) - xyz_QvapSat(i,j,k-1) ) + xyz_DPressDz(i,j,k) * (xyz_Qvap(i,j,k) - xyz_QvapSat(i,j,k) )
DQvapSatDTempUpper = EL * xyz_QvapSat(i,j,k) / ( RVap * xyz_Temp(i,j,k) * xyz_Temp(i,j,k) )
DQvapSatDTempUnder = EL * xyz_QvapSat(i,j,k-1) / ( RVap * xyz_Temp(i,j,k-1) * xyz_Temp(i,j,k-1) )
DHDTempUpper = 1. + EL/Cp * DQvapSatDTempUpper
DHDTempUnder = 1. + EL/Cp * DQvapSatDTempUnder
DelTempSat = TempSat + ( 1. - xyz_DPFact(i,j,k) / DHDTempUnder ) * EL/Cp * DelQvap / xyz_DPressDz(i,j,k-1)
! 温度の調節
DelTempUpper = DelTempSat / ( ( 1. + xyz_DDPressDDPress(i,j,k) ) * DHDTempUpper + xyz_DPFact(i,j,k) * ( 1. - xyz_DDPressDDPress(i,j,k) * DHDTempUpper / DHDTempUnder ) )
DelTempUnder = - DHDTempUpper / DHDTempUnder * xyz_DDPressDDPress(i,j,k) * DelTempUpper + EL/Cp * DelQvap / ( xyz_DPressDz(i,j,k-1) * DHDTempUnder )
xyz_Temp(i,j,k) = xyz_Temp(i,j,k) + DelTempUpper
xyz_Temp(i,j,k-1) = xyz_Temp(i,j,k-1) + DelTempUnder
! 比湿の調節
xyz_Qvap(i,j,k) = xyz_QvapSat(i,j,k) + DQvapSatDTempUpper * DelTempUpper
xyz_Qvap(i,j,k-1) = xyz_QvapSat(i,j,k-1) + DQvapSatDTempUnder * DelTempUnder
xyz_QvapSat(i,j,k) = xyz_Qvap(i,j,k)
xyz_QvapSat(i,j,k-1) = xyz_Qvap(i,j,k-1)
! 調節したか否か
xy_Adjust(i,j) = 1.
end if
end if
end if
end do
end do
end do
Adjust = 0.
do i = 1, im
do j = 1, jm
xy_Adjust_b(i,j) = xy_Adjust(i,j)
Adjust = Adjust + xy_Adjust(i,j)
end do
end do
if ( Adjust .LT. 1. ) exit
end do
!----- 比湿変化率, 温度変化率, 降水量の算出 -----
do k = 1, km
do i = 1, im
do j = 1, jm
! 比湿変化率
xyz_DCumulusQvapDt(i,j,k) = xyz_DCumulusQvapDt(i,j,k) + ( xyz_Qvap(i,j,k) - xyz_Qvap_b(i,j,k) ) / DelTimePhy
! 温度変化率
xyz_DCumulusTempDt(i,j,k) = xyz_DCumulusTempDt(i,j,k) + ( xyz_Temp(i,j,k) - xyz_Temp_b(i,j,k) ) / DelTimePhy
! 降水量
xy_CumulusRain(i,j) = xy_CumulusRain(i,j) + ( xyz_Qvap_b(i,j,k) - xyz_Qvap(i,j,k) ) * xyz_DPressDz(i,j,k) / Grav * EL / DelTimePhy
end do
end do
end do
!----------------------------------------------------------------
! 終了処理
!----------------------------------------------------------------
call EndSub(subname)
end subroutine physics_cumulus_adjust