subroutine physics_verdiff_main( xyr_VelLonFlux , xyr_VelLatFlux , xyr_TempFlux , xyr_QvapFlux , xyzo_VelMatrix , xyzo_TempMatrix , xyzo_QvapMatrix , xyz_VelLon , xyz_VelLat , xyz_Temp , xyr_Temp , xyz_Qvap , xyz_Press , xyr_Press , xyz_GeoPot , xyr_GeoPot ) ! (in) 高度 (半整数)
!==== Dependency
use type_mod, only: REKIND, DBKIND, INTKIND, TOKEN, STRING
use grid_3d_mod, only: im, jm, km
use constants_mod, only: RAir, Cp, Grav, EL
use physics_verdiff_coeff_mod, only: physics_verdiff_coeff
use dc_trace, only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump
implicit none
!==== Output
!
real(DBKIND), intent(out) :: xyr_VelLonFlux(im,jm,km+1) , xyr_VelLatFlux(im,jm,km+1) , xyr_TempFlux(im,jm,km+1) , xyr_QvapFlux(im,jm,km+1) ! (out) 比湿フラックス
!==== Input
!
real(DBKIND), intent(in) :: xyz_VelLon(im,jm,km) , xyz_VelLat(im,jm,km) , xyz_Temp(im,jm,km) , xyr_Temp(im,jm,km+1) , xyz_Qvap(im,jm,km) , xyz_Press(im,jm,km) , xyr_Press(im,jm,km+1) , xyz_GeoPot(im,jm,km) , xyr_GeoPot(im,jm,km+1) ! (in) 高度 (半整数)
!==== In/Out
!
real(DBKIND), intent(inout) :: xyzo_VelMatrix(im,jm,km,-1:1) , xyzo_TempMatrix(im,jm,0:km,-1:1) , xyzo_QvapMatrix(im,jm,km,-1:1) ! (inout) 比湿陰解行列
!----- 作業用内部変数 -----
character(STRING), parameter:: subname = "physics_verdiff_main"
! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
integer(INTKIND) :: i,j, k
real(DBKIND), parameter :: RefPress = 100000.0D0 ! 参照気圧
real(DBKIND), parameter :: BasePotTemp = 300.0D0 ! 基本温位
real(DBKIND), parameter :: SquareVelMin = 0.1D0 ! 風二乗差最小値
real(DBKIND), parameter :: BulkRiNumMin = - 100.0D0 ! バルクＲｉ数最小値
real(DBKIND) :: xyr_DVelDz(im,jm,km+1) , xyr_BulkRiNum(im,jm,km+1) , xyr_TempTransCoeff(im,jm,km+1) , xyr_QvapTransCoeff(im,jm,km+1) , xyr_VelTransCoeff(im,jm,km+1) , xyr_TempDiffuCoeff(im,jm,km+1) , xyr_QvapDiffuCoeff(im,jm,km+1) , xyr_VelDiffuCoeff(im,jm,km+1) , xyz_Exner(im,jm,km) , xyr_Exner(im,jm,km+1) ! Exner関数 (半整数)
continue
!----------------------------------------------------------------
! 開始処理
!----------------------------------------------------------------
call BeginSub(subname)
!----------------------------------------------------------------
! 鉛直拡散計算
!----------------------------------------------------------------
! ---- 1. 定数計算 ----
xyz_Exner = ( xyz_Press / RefPress ) ** (RAir/Cp)
xyr_Exner = ( xyr_Press / RefPress ) ** (RAir/Cp)
! ---- 2. バルクＲｉ数 ----
xyr_DVelDz = 0.0d0
xyr_BulkRiNum = 0.0d0
do k = 2, km
do i = 1, im
do j = 1, jm
xyr_DVelDz(i,j,k) = SQRT( MAX( SquareVelMin , ( xyz_VelLon(i,j,k) - xyz_VelLon(i,j,k-1) )**2 + ( xyz_VelLat(i,j,k) - xyz_VelLat(i,j,k-1) )**2 ) ) / ( xyz_GeoPot(i,j,k) - xyz_GeoPot(i,j,k-1) )
xyr_BulkRiNum(i,j,k) = Grav / BasePotTemp * ( xyz_Temp(i,j,k) / xyz_Exner(i,j,k) - xyz_Temp(i,j,k-1) / xyz_Exner(i,j,k-1) ) / ( xyz_GeoPot(i,j,k) - xyz_GeoPot(i,j,k-1) ) / xyr_DVelDz(i,j,k)**2
xyr_BulkRiNum(i,j,k) = MAX( xyr_BulkRiNum(i,j,k) , BulkRiNumMin )
end do
end do
end do
! ---- 3. 拡散係数 ----
call physics_verdiff_coeff( xyr_VelDiffuCoeff , xyr_TempDiffuCoeff, xyr_QvapDiffuCoeff, xyr_BulkRiNum , xyr_DVelDz , xyr_GeoPot ) ! (in)
! ---- *. 浅い積雲対流 ----
! ---- *. 拡散係数の出力 ----
! CALL HISTIN ( DFM , 'DFM' )
! CALL HISTIN ( DFH , 'DFH' )
! CALL HISTIN ( DFE , 'DFE' )
! ---- 5. 輸送係数 ----
xyr_VelTransCoeff = 0.0d0
xyr_TempTransCoeff = 0.0d0
xyr_QvapTransCoeff = 0.0d0
do k = 2, km
xyr_VelTransCoeff(:,:,k) = xyr_VelDiffuCoeff(:,:,k) * xyr_Press(:,:,k) / RAir / xyr_Temp(:,:,k) / ( xyz_GeoPot(:,:,k) - xyz_GeoPot(:,:,k-1) )
xyr_TempTransCoeff(:,:,k) = xyr_TempDiffuCoeff(:,:,k) * xyr_Press(:,:,k) / RAir / xyr_Temp(:,:,k) / ( xyz_GeoPot(:,:,k) - xyz_GeoPot(:,:,k-1) )
xyr_QvapTransCoeff(:,:,k) = xyr_QvapDiffuCoeff(:,:,k) * xyr_Press(:,:,k) / RAir / xyr_Temp(:,:,k) / ( xyz_GeoPot(:,:,k) - xyz_GeoPot(:,:,k-1) )
end do
! ---- 5. フラックス----
! k =0, km+1 も 0 でいいのかな.
xyr_VelLonFlux = 0.0d0
xyr_VelLatFlux = 0.0d0
xyr_TempFlux = 0.0d0
xyr_QvapFlux = 0.0d0
do k = 2, km
xyr_VelLonFlux(:,:,k) = xyr_VelLonFlux(:,:,k) + xyr_VelTransCoeff(:,:,k) * ( xyz_VelLon(:,:,k-1) - xyz_VelLon(:,:,k) )
xyr_VelLatFlux(:,:,k) = xyr_VelLatFlux(:,:,k) + xyr_VelTransCoeff(:,:,k) * ( xyz_VelLat(:,:,k-1) - xyz_VelLat(:,:,k) )
xyr_TempFlux(:,:,k) = xyr_TempFlux(:,:,k) + Cp * xyr_TempTransCoeff(:,:,k) * xyr_Exner(:,:,k) * ( xyz_Temp(:,:,k-1) / xyz_Exner(:,:,k-1) - xyz_Temp(:,:,k) / xyz_Exner(:,:,k) )
xyr_QvapFlux(:,:,k) = xyr_QvapFlux(:,:,k) + EL * xyr_QvapTransCoeff(:,:,k) * ( xyz_Qvap(:,:,k-1) - xyz_Qvap(:,:,k) )
end do
! ---- 5. 陰解用行列 ----
do k = 2, km
xyzo_VelMatrix(:,:,k,0) = xyzo_VelMatrix(:,:,k,0) + xyr_VelTransCoeff(:,:,k)
xyzo_VelMatrix(:,:,k,-1) = - xyr_VelTransCoeff(:,:,k)
xyzo_TempMatrix(:,:,k,0) = xyzo_TempMatrix(:,:,k,0) + Cp * xyr_TempTransCoeff(:,:,k) * xyr_Exner(:,:,k) / xyz_Exner(:,:,k)
xyzo_TempMatrix(:,:,k,-1) = - Cp * xyr_TempTransCoeff(:,:,k) * xyr_Exner(:,:,k) / xyz_Exner(:,:,k-1)
xyzo_QvapMatrix(:,:,k,0) = xyzo_QvapMatrix(:,:,k,0) + Cp * xyr_QvapTransCoeff(:,:,k)
xyzo_QvapMatrix(:,:,k,-1) = - Cp * xyr_QvapTransCoeff(:,:,k)
end do
do k = 1, km-1
xyzo_VelMatrix(:,:,k,0) = xyzo_VelMatrix(:,:,k,0) + xyr_VelTransCoeff(:,:,k+1)
xyzo_VelMatrix(:,:,k,1) = - xyr_VelTransCoeff(:,:,k+1)
xyzo_TempMatrix(:,:,k,0) = xyzo_TempMatrix(:,:,k,0) + Cp * xyr_TempTransCoeff(:,:,k+1) * xyr_Exner(:,:,k+1) / xyz_Exner(:,:,k)
xyzo_TempMatrix(:,:,k,1) = - Cp * xyr_TempTransCoeff(:,:,k+1) * xyr_Exner(:,:,k+1) / xyz_Exner(:,:,k+1)
xyzo_QvapMatrix(:,:,k,0) = xyzo_QvapMatrix(:,:,k,0) + Cp * xyr_QvapTransCoeff(:,:,k+1)
xyzo_QvapMatrix(:,:,k,1) = - Cp * xyr_QvapTransCoeff(:,:,k+1)
end do
!----------------------------------------------------------------
! 終了処理
!----------------------------------------------------------------
call EndSub(subname)
end subroutine physics_verdiff_main