subroutine physics_verdiff_coeff( xyr_VelDiffuCoeff , xyr_TempDiffuCoeff, xyr_QvapDiffuCoeff, xyr_BulkRiNum , xyr_DVelDz , 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: FKarm
use dc_trace, only: SetDebug, BeginSub, EndSub, DbgMessage, DataDump
implicit none
!==== Output
!
real(DBKIND), intent(out) :: xyr_TempDiffuCoeff(im,jm,km+1) , xyr_QvapDiffuCoeff(im,jm,km+1) , xyr_VelDiffuCoeff(im,jm,km+1) ! (out) 拡散係数：運動量
!==== Input
!
real(DBKIND), intent(in) :: xyr_DVelDz(im,jm,km+1) , xyr_BulkRiNum(im,jm,km+1) , xyr_GeoPot(im,jm,km+1) ! (in) 高度 (半整数)
!----- 作業用内部変数 -----
character(STRING), parameter:: subname = "physics_verdiff_coeff"
! do ループ用作業変数 (東西 i*、南北 j*、鉛直 k*、波数 l*用)
integer(INTKIND) :: i, j, k
real(DBKIND), parameter :: MixLengthMax = 300.0d0 ! 最大混合距離
real(DBKIND), parameter :: TildeShMin = 0.0d0 ! tilde(Sh)最小値
real(DBKIND), parameter :: TildeSmMin = 0.0d0 ! tilde(Sm)最小値
real(DBKIND), parameter :: VelDiffuCoeffMin = 0.1d0 ! u 拡散係数最小値
real(DBKIND), parameter :: TempDiffuCoeffMin = 0.1d0 ! T 拡散係数最小値
real(DBKIND), parameter :: QvapDiffuCoeffMin = 0.1d0 ! q 拡散係数最小値
real(DBKIND), parameter :: VelDiffuCoeffMax = 10000.0d0 ! u 拡散係数最大値
real(DBKIND), parameter :: TempDiffuCoeffMax = 10000.0d0 ! T 拡散係数最大値
real(DBKIND), parameter :: QvapDiffuCoeffMax = 10000.0d0 ! q 拡散係数最大値
real(DBKIND), parameter :: ParamA1 = 0.92 ! Mellor Yamada Lev2 定数
real(DBKIND), parameter :: ParamB1 = 16.6 ! Mellor Yamada Lev2 定数
real(DBKIND), parameter :: ParamA2 = 0.74 ! Mellor Yamada Lev2 定数
real(DBKIND), parameter :: ParamB2 = 10.1 ! Mellor Yamada Lev2 定数
real(DBKIND), parameter :: ParamC1 = 0.08 ! Mellor Yamada Lev2 定数
real(DBKIND) :: xyr_FluxRiNum(im,jm,km+1) , xyr_TildeSh(im,jm,km+1) , xyr_TildeSm(im,jm,km+1) , xyr_MixLength(im,jm,km+1) ! 混合距離
real(DBKIND) :: Alpha1, Alpha2, Beta1, Beta2, Beta3, Beta4, Gamma1, Gamma2, CriticalFluxRiNum
continue
!----------------------------------------------------------------
! 開始処理
!----------------------------------------------------------------
call BeginSub(subname)
!----------------------------------------------------------------
! 鉛直拡散計算
!----------------------------------------------------------------
! ---- 1. 定数計算 ----
Gamma1 = ( 1./ 3.) - ( 2.* ParamA1 / ParamB1 )
Gamma2 = ( ParamB2 / ParamB1 ) + ( 6.* ParamA1 / ParamB1 )
Alpha1 = 3. * ParamA2 * Gamma1
Alpha2 = 3. * ParamA2 * ( Gamma1 + Gamma1 )
Beta1 = ParamA1 * ParamB1 * ( Gamma1 - ParamC1 )
Beta2 = ParamA1 * ( ParamB1 * ( Gamma1 - ParamC1 ) + 6.*ParamA1 + 3.*ParamA2 )
Beta3 = ParamA2 * ParamB1 * Gamma1
Beta4 = ParamA2 * ( ParamB1 * ( Gamma1 + Gamma2 ) - 3.*ParamA1 )
CriticalFluxRiNum = Gamma1 / ( Gamma1 + Gamma2 )
! ---- 2. フラックスＲｉ数 ----
xyr_FluxRiNum = ( Beta1 + Beta4 * xyr_BulkRiNum - SQRT( ( Beta1 + Beta4 * xyr_BulkRiNum )**2 - 4. * Beta2 * Beta3 * xyr_BulkRiNum ) ) / ( 2. * Beta2 )
! ---- 3. tilde(Sm), tilde(Sh) ----
xyr_TildeSh = 0.0d0
xyr_TildeSm = 0.0d0
do k = 1, km
do i = 1, im
do j = 1, jm
if ( xyr_FluxRiNum(i,j,k) .LT. CriticalFluxRiNum ) then
xyr_TildeSh(i,j,k) = ( Alpha1 - Alpha2 * xyr_FluxRiNum(i,j,k) ) / ( 1. - xyr_FluxRiNum(i,j,k) )
xyr_TildeSm(i,j,k) = ( Beta1 - Beta2 * xyr_FluxRiNum(i,j,k) ) / ( Beta3 - Beta4 * xyr_FluxRiNum(i,j,k) ) * xyr_TildeSh(i,j,k)
xyr_TildeSh(i,j,k) = MAX ( xyr_TildeSh(i,j,k) , TildeShMin )
xyr_TildeSm(i,j,k) = MAX ( xyr_TildeSm(i,j,k) , TildeSmMin )
else
xyr_TildeSh(i,j,k) = TildeShMin
xyr_TildeSm(i,j,k) = TildeSmMin
end if
end do
end do
end do
! ---- 4. 混合距離 ----
xyr_MixLength = FKarm * xyr_GeoPot / (1. + FKarm * xyr_GeoPot / MixLengthMax )
! ---- 5. 拡散係数 ----
xyr_VelDiffuCoeff = xyr_MixLength ** 2 * xyr_DVelDz * SQRT ( ParamB1 * ( 1. - xyr_FluxRiNum ) * xyr_TildeSm ) * xyr_TildeSm
xyr_TempDiffuCoeff = xyr_MixLength ** 2 * xyr_DVelDz * SQRT ( ParamB1 * ( 1. - xyr_FluxRiNum ) * xyr_TildeSm ) * xyr_TildeSh
xyr_QvapDiffuCoeff = xyr_TempDiffuCoeff
do k = 1, km
do i = 1, im
do j = 1, jm
xyr_VelDiffuCoeff(i,j,k) = MAX( MIN( xyr_VelDiffuCoeff(i,j,k) , VelDiffuCoeffMax ) , VelDiffuCoeffMin )
xyr_TempDiffuCoeff(i,j,k) = MAX( MIN( xyr_TempDiffuCoeff(i,j,k) , TempDiffuCoeffMax ) , TempDiffuCoeffMin )
xyr_QvapDiffuCoeff(i,j,k) = MAX( MIN( xyr_QvapDiffuCoeff(i,j,k) , QvapDiffuCoeffMax ) , QvapDiffuCoeffMin )
end do
end do
end do
xyr_VelDiffuCoeff(:,:,1) = 0.0d0
xyr_TempDiffuCoeff(:,:,1) = 0.0d0
xyr_QvapDiffuCoeff(:,:,1) = 0.0d0
xyr_VelDiffuCoeff(:,:,km+1) = 0.0d0
xyr_TempDiffuCoeff(:,:,km+1) = 0.0d0
xyr_QvapDiffuCoeff(:,:,km+1) = 0.0d0
!----------------------------------------------------------------
! 終了処理
!----------------------------------------------------------------
call EndSub(subname)
end subroutine physics_verdiff_coeff