Subroutine : |
|
pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(in)
|
xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(in)
|
xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(in)
|
xyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(in)
|
xyz_ExnerBl(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(in)
|
xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax, ncmax) : | real(DP),intent(in)
|
xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(in)
|
xyz_KhBl(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(in)
|
xyz_CDensBl(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(in)
|
pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(inout)
|
xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(inout)
|
xyr_DVelZDt(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(inout)
|
xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(inout)
|
xyz_DExnerDt(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(inout)
|
xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax) : | real(DP),intent(inout)
|
xyz_DKmDt(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(inout)
|
xyz_DCDensDt(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(inout)
|
xyz_KmAl(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(out)
|
xyz_KhAl(imin:imax,jmin:jmax,kmin:kmax) : | real(DP),intent(out)
|
subroutine turbulence_KW1978_forcing( pyz_VelXBl, xqz_VelYBl, xyr_VelZBl, xyz_PTempBl, xyz_ExnerBl, xyzf_QMixBl, xyz_KmBl, xyz_KhBl, xyz_CDensBl, pyz_DVelXDt, xqz_DVelYDt, xyr_DVelZDt, xyz_DPTempDt,xyz_DExnerDt, xyzf_DQMixDt, xyz_DKmDt, xyz_DCDensDt, xyz_KmAl, xyz_KhAl )
implicit none
real(DP),intent(in) :: pyz_VelXBl(imin:imax,jmin:jmax,kmin:kmax)
!水平速度
real(DP),intent(in) :: xqz_VelYBl(imin:imax,jmin:jmax,kmin:kmax)
!水平速度
real(DP),intent(in) :: xyr_VelZBl(imin:imax,jmin:jmax,kmin:kmax)
!鉛直速度
real(DP),intent(in) :: xyz_PTempBl(imin:imax,jmin:jmax,kmin:kmax)
!温位
real(DP),intent(in) :: xyz_ExnerBl(imin:imax,jmin:jmax,kmin:kmax)
!無次元圧力
real(DP),intent(in) :: xyzf_QMixBl(imin:imax,jmin:jmax,kmin:kmax, ncmax)
!凝縮成分の混合比
real(DP),intent(in) :: xyz_KmBl(imin:imax,jmin:jmax,kmin:kmax)
!乱流拡散係数
real(DP),intent(in) :: xyz_KhBl(imin:imax,jmin:jmax,kmin:kmax)
!乱流拡散係数
real(DP),intent(in) :: xyz_CDensBl(imin:imax,jmin:jmax,kmin:kmax)
real(DP),intent(inout):: pyz_DVelXDt(imin:imax,jmin:jmax,kmin:kmax)
!スカラー量の水平乱流拡散
real(DP),intent(inout):: xqz_DVelYDt(imin:imax,jmin:jmax,kmin:kmax)
!スカラー量の水平乱流拡散
real(DP),intent(inout):: xyr_DVelZDt(imin:imax,jmin:jmax,kmin:kmax)
!スカラー量の水平乱流拡散
real(DP),intent(inout):: xyz_DPTempDt(imin:imax,jmin:jmax,kmin:kmax)
real(DP),intent(inout):: xyz_DExnerDt(imin:imax,jmin:jmax,kmin:kmax)
real(DP),intent(inout):: xyzf_DQMixDt(imin:imax,jmin:jmax,kmin:kmax, ncmax)
real(DP),intent(inout):: xyz_DKmDt(imin:imax,jmin:jmax,kmin:kmax)
real(DP),intent(inout) :: xyz_DCDensDt(imin:imax,jmin:jmax,kmin:kmax)
real(DP),intent(out):: xyz_KmAl(imin:imax,jmin:jmax,kmin:kmax)
!乱流拡散係数
real(DP),intent(out):: xyz_KhAl(imin:imax,jmin:jmax,kmin:kmax)
!乱流拡散係数
real(DP) :: xyz_Buoy(imin:imax,jmin:jmax,kmin:kmax)
!渦粘性係数の
real(DP) :: xyz_Shear(imin:imax,jmin:jmax,kmin:kmax)
!渦粘性係数の
real(DP) :: xyz_Disp(imin:imax,jmin:jmax,kmin:kmax)
!乱流エネルギーの消散
real(DP) :: xyz_DispPI(imin:imax,jmin:jmax,kmin:kmax)
!乱流エネルギーの消散
real(DP) :: xyz_DispHeat(imin:imax,jmin:jmax,kmin:kmax)
!乱流エネルギーの消散
real(DP) :: xyz_Turb(imin:imax,jmin:jmax,kmin:kmax)
!
real(DP) :: xyzf_Turb(imin:imax,jmin:jmax,kmin:kmax, ncmax)
!スカラー量の水平乱流拡散
real(DP) :: pyz_Turb(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xqz_Turb(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyr_Turb(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: pyz_DVelXDt0(imin:imax,jmin:jmax,kmin:kmax)
!スカラー量の水平乱流拡散
real(DP) :: xqz_DVelYDt0(imin:imax,jmin:jmax,kmin:kmax)
!スカラー量の水平乱流拡散
real(DP) :: xyr_DVelZDt0(imin:imax,jmin:jmax,kmin:kmax)
!スカラー量の水平乱流拡散
real(DP) :: xyz_DPTempDt0(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyz_DExnerDt0(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyzf_DQMixDt0(imin:imax,jmin:jmax,kmin:kmax, ncmax)
real(DP) :: xyz_DKmDt0(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyz_DCDensDt0(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyz_PTempBlAll(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyzf_QMixBlAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
integer :: f
!----------------------------------
! 初期化
!
xyz_PTempBlAll = xyz_PTempBl + xyz_PTempBZ
xyzf_QMixBlAll = xyzf_QMixBl + xyzf_QMixBZ
pyz_DVelXDt0 = pyz_DVelXDt
xqz_DVelYDt0 = xqz_DVelYDt
xyr_DVelZDt0 = xyr_DVelZDt
xyz_DPTempDt0 = xyz_DPTempDt
xyz_DExnerDt0 = xyz_DExnerDt
xyzf_DQMixDt0 = xyzf_DQMixDt
xyz_DKmDt0 = xyz_DKmDt
xyz_DCDensDt0 = xyz_DCDensDt
!----------------------------------
! 拡散係数の時間発展 (エネルギー方程式を Km の式に変形したもの)
!
! Buoyancy term
!
xyz_Buoy = xyz_BuoyMoistKm(xyz_PTempBl, xyz_ExnerBl, xyzf_QMixBl)
! xyz_Buoy = &
! & - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) &
! & * xyz_avr_xyr( xyr_dz_xyz( xyz_PTempBlAll ) ) &
! & / ( 2.0d0 * xyz_PTempBZ )
xyz_Shear = ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) * ( ( xyz_dx_pyz( pyz_VelXBl ) ) ** 2.0d0 + ( xyz_dy_xqz( xqz_VelYBl ) ) ** 2.0d0 + ( xyz_dz_xyr( xyr_VelZBl ) ) ** 2.0d0 + 5.0d-1 * ( ( xyz_avr_pyr( pyr_dz_pyz( pyz_VelXBl ) ) + xyz_avr_pyr( pyr_dx_xyr( xyr_VelZBl ) ) ) ** 2.0d0 + ( xyz_avr_xqr( xqr_dy_xyr( xyr_VelZBl ) ) + xyz_avr_xqr( xqr_dz_xqz( xqz_VelYBl ) ) ) ** 2.0d0 + ( xyz_avr_pqz( pqz_dx_xqz( xqz_VelYBl ) ) + xyz_avr_pqz( pqz_dy_pyz( pyz_VelXBl ) ) ) ** 2.0d0 ) ) - xyz_KmBl * ( xyz_dx_pyz( pyz_VelXBl ) + xyz_dy_xqz( xqz_VelYBl ) + xyz_dz_xyr( xyr_VelZBl ) ) / 3.0d0
xyz_Disp = - (xyz_KmBl ** 2.0d0) * 5.0d-1 / (MixLen ** 2.0d0)
! tendency
!
xyz_DKmDt = xyz_DKmDt0 + xyz_Buoy + xyz_Shear + xyz_Disp
! 時間積分
!
xyz_KmAl = xyz_KmBl + 2.0d0 * DelTimeLong * xyz_DKmDt
! 上限値の設定
!
xyz_KmAl = max( 0.0d0, min( xyz_KmAl, KmMax ) )
! Kh の設定
!
xyz_KhAl = 3.0d0 * xyz_KmAl
!--------------------------------
! 雲密度の tendency
!
xyz_Turb = xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyz_CDensBl ) ) + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyz_CDensBl ) ) + xyz_dz_xyr( xyr_avr_xyz( xyz_KhBl ) * xyr_dz_xyz( xyz_CDensBl ) )
xyz_DCDensDt = xyz_DCDensDt0 + xyz_Turb
! output
!
call HistoryAutoPut(TimeN, 'CDensTurb', xyz_Turb(1:nx,1:ny,1:nz))
!--------------------------------
! 温位の tendency
!
xyz_Turb = xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyz_PTempBlAll ) ) + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyz_PTempBlAll ) ) + xyz_dz_xyr( xyr_avr_xyz( xyz_KhBl ) * xyr_dz_xyz( xyz_PTempBlAll ) )
xyz_DispHeat = (xyz_KmBl ** 3.0d0) / (xyz_ExnerBZ * CpDry * (Cm ** 2.0d0) * (MixLen ** 4.0d0))
xyz_DPTempDt = xyz_DPTempDt0 + xyz_Turb + xyz_DispHeat
call HistoryAutoPut(TimeN, 'PTempDisp', xyz_DispHeat(1:nx, 1:ny, 1:nz))
call HistoryAutoPut(TimeN, 'PTempTurb', xyz_Turb(1:nx, 1:ny, 1:nz))
!--------------------------------
! 混合比の tendency
!
do f = 1, ncmax
xyzf_Turb(:,:,:,f) = xyz_dx_pyz( pyz_avr_xyz( xyz_KhBl ) * pyz_dx_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) + xyz_dy_xqz( xqz_avr_xyz( xyz_KhBl ) * xqz_dy_xyz( xyzf_QMixBlAll(:,:,:,f) ) ) + xyz_dz_xyr( xyr_avr_xyz( xyz_KhBl ) * xyr_dz_xyz( xyzf_QMixBlAll(:,:,:,f) ) )
call HistoryAutoPut(TimeN, trim(SpcWetSymbol(f))//'_Turb', xyzf_Turb(1:nx,1:ny,1:nz,f))
end do
xyzf_DQMixDt = xyzf_DQMixDt0 + xyzf_Turb
!--------------------------------
! 各速度成分の tendency
!
pyz_Turb = 2.0d0 * pyz_dx_xyz( xyz_KmBl * xyz_dx_pyz( pyz_VelXBl ) ) + pyz_dy_pqz( pqz_avr_xyz( xyz_KmBl ) * pqz_dx_xqz( xqz_VelYBl ) + pqz_avr_xyz( xyz_KmBl ) * pqz_dy_pyz( pyz_VelXBl ) ) + pyz_dz_pyr( pyr_avr_xyz( xyz_KmBl ) * pyr_dx_xyr( xyr_VelZBl ) + pyr_avr_xyz( xyz_KmBl ) * pyr_dz_pyz( pyz_VelXBl ) ) - 2.0d0 * pyz_dx_xyz( ( xyz_KmBl ** 2.0d0 ) ) / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )
pyz_DVelXDt = pyz_DVelXDt0 + pyz_Turb
call HistoryAutoPut(TimeN, 'VelXTurb', pyz_Turb(1:nx, 1:ny, 1:nz))
if (jmin == jmax) then
xqz_DVelYDt = xqz_DVelYDt0
else
xqz_Turb = 2.0d0 * xqz_dy_xyz( xyz_KmBl * xyz_dy_xqz( xqz_VelYBl ) ) + xqz_dx_pqz( pqz_avr_xyz( xyz_KmBl ) * pqz_dy_pyz( pyz_VelXBl ) + pqz_avr_xyz( xyz_KmBl ) * pqz_dx_xqz( xqz_VelYBl ) ) + xqz_dz_xqr( xqr_avr_xyz( xyz_KmBl ) * xqr_dy_xyr( xyr_VelZBl ) + xqr_avr_xyz( xyz_KmBl ) * xqr_dz_xqz( xqz_VelYBl ) ) - 2.0d0 * xqz_dy_xyz( ( xyz_KmBl ** 2.0d0 ) ) / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )
xqz_DVelYDt = xqz_DVelYDt0 + xqz_Turb
end if
call HistoryAutoPut(TimeN, 'VelYTurb', xqz_Turb(1:nx, 1:ny, 1:nz))
xyr_Turb = + 2.0d0 * xyr_dz_xyz( xyz_KmBl * xyz_dz_xyr( xyr_VelZBl ) ) + xyr_dx_pyr( pyr_avr_xyz( xyz_KmBl ) * pyr_dz_pyz( pyz_VelXBl ) + pyr_avr_xyz( xyz_KmBl ) * pyr_dx_xyr( xyr_VelZBl ) ) + xyr_dy_xqr( xqr_avr_xyz( xyz_KmBl ) * xqr_dz_xqz( xqz_VelYBl ) + xqr_avr_xyz( xyz_KmBl ) * xqr_dy_xyr( xyr_VelZBl ) ) - 2.0d0 * xyr_dz_xyz( ( xyz_KmBl ** 2.0d0 ) ) / ( 3.0d0 * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) )
xyr_DVelZDt = xyr_DVelZDt0 + xyr_Turb
call HistoryAutoPut(TimeN, 'VelZTurb', xyr_Turb(1:nx, 1:ny, 1:nz))
!--------------------
! Exner function
!
xyz_DispPI = (xyz_VelSoundBZ ** 2.0d0) * xyz_DispHeat / (CpDry * xyz_VPTempBZ ** 2.0d0)
xyz_DExnerDt = xyz_DExnerDt0 + xyz_DispPI
call HistoryAutoPut(TimeN, 'ExnerDisp', xyz_DispPI(1:nx, 1:ny, 1:nz))
end subroutine Turbulence_KW1978_forcing