subroutine Cloudphys_K1969_forcing(xyz_ExnerNl, xyz_DExnerDt, xyz_PTempAl, xyzf_QMixAl)
implicit none
real(DP), intent(in) :: xyz_ExnerNl(imin:imax, jmin:jmax, kmin:kmax)
real(DP), intent(inout) :: xyz_DExnerDt(imin:imax, jmin:jmax, kmin:kmax)
real(DP), intent(inout) :: xyz_PTempAl(imin:imax, jmin:jmax, kmin:kmax)
real(DP), intent(inout) :: xyzf_QMixAl(imin:imax, jmin:jmax, kmin:kmax, ncmax)
real(DP) :: xyz_PTempOrig(imin:imax, jmin:jmax, kmin:kmax)
real(DP) :: xyz_PTempWork(imin:imax, jmin:jmax, kmin:kmax)
real(DP) :: xyz_DelPTemp(imin:imax, jmin:jmax, kmin:kmax)
real(DP) :: xyz_PTempCond(imin:imax, jmin:jmax, kmin:kmax)
real(DP) :: xyzf_QMixOrig(imin:imax, jmin:jmax, kmin:kmax, ncmax)
real(DP) :: xyzf_QMixWork(imin:imax, jmin:jmax, kmin:kmax, ncmax)
real(DP) :: xyzf_DelQMix(imin:imax, jmin:jmax, kmin:kmax, ncmax)
real(DP) :: xyzf_QMixCond(imin:imax, jmin:jmax, kmin:kmax, ncmax)
real(DP) :: DelTime
integer :: l, s
integer :: iG, iC, iR
real(DP) :: xyzf_Cloud2Rain(imin:imax,jmin:jmax,kmin:kmax, ncmax)
!雲から雨への変換量
real(DP) :: xyz_AutoConv(imin:imax,jmin:jmax,kmin:kmax)
!飽和混合比
real(DP) :: xyz_Collect(imin:imax,jmin:jmax,kmin:kmax)
!規格化された潜熱
real(DP) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
!混合比の擾乱成分 + 平均成分
real(DP) :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax)
!温度の擾乱成分 + 平均成分
real(DP) :: xyz_PressAll(imin:imax,jmin:jmax,kmin:kmax)
!全圧
real(DP) :: xyz_ExnerAll(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyz_NonSaturate(imin:imax,jmin:jmax,kmin:kmax)
!未飽和度(飽和混合比と蒸気の混合比の差)
real(DP) :: xyzf_Rain2Gas(imin:imax,jmin:jmax,kmin:kmax, ncmax)
real(DP) :: xyzf_Rain2GasNH4SH(imin:imax,jmin:jmax,kmin:kmax, ncmax)
real(DP) :: xyzf_DelPTemp(imin:imax,jmin:jmax,kmin:kmax, ncmax)
real(DP) :: xyz_DelPTempNH4SH(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyz_PressDry(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyz_ExnerCond(imin:imax,jmin:jmax,kmin:kmax)
!-----------------------------------------
! 時間刻み幅. Leap-frog なので, 2 \del t
!
DelTime = 2.0d0 * DelTimeLong
!------------------------------------------
! 初期値を保管 Store Initial Value
!
xyz_PTempOrig = xyz_PTempAl
xyzf_QMixOrig = xyzf_QMixAl
! 全エクスナー関数・全圧を計算. サブルーチン内では変化しない.
!
xyz_ExnerAll = xyz_ExnerNl + xyz_ExnerBZ
xyz_PressAll = PressBasis * (xyz_ExnerAll ** (CpDry / GasRDry))
!------------------------------------------
! 暖かい雨のパラメタリゼーション.
! * 雲<-->雨 の変換を行う.
!
! Warm rain parameterization.
! * Conversion from cloud to rain.
!これまでの値を作業配列に保管
! Previous values are stored to work area.
!
xyzf_QMixWork = xyzf_QMixAl
!雨への変化量を計算
! Conversion values are calculated.
!
xyzf_QMixAll = max( 1.0d-60, xyzf_QMixAl + xyzf_QMixBZ )
xyzf_Cloud2Rain = 0.0d0
do s = 1, CloudNum
! 値を保管
!
iC = IdxC(s)
iR = IdxR(s)
!併合成長
!
xyz_AutoConv = DelTime / AutoConvTime * max( 1.0d-60, ( xyzf_QMixAll(:,:,:,iC) - QMixCr) )
!衝突合体成長
!
xyz_Collect = DelTime * 2.2d0 * FactorJ * xyzf_QMixAll(:,:,:,iC) * (xyzf_QMixAll(:,:,:,iR) * xyz_DensBZ) ** 0.875d0
!雲の変換量: 併合成長と合体衝突の和
! 元々の変化量を上限値として設定する. 負の値となる.
!
xyzf_Cloud2Rain(:,:,:,iC) = - min( xyzf_QMixAll(:,:,:,iC), ( xyz_AutoConv + xyz_Collect ) )
!雨の変換量. 符号は雲の変換量とは反対.
xyzf_Cloud2Rain(:,:,:,iR) = - xyzf_Cloud2Rain(:,:,:,iC)
end do
! 変化量を足し込む
!
xyzf_QMixAl = xyzf_QMixWork + xyzf_Cloud2Rain
! Set Margin
!
! call SetMargin_xyzf(xyzf_QMixAl)
!-------------------------------------------
! 湿潤飽和調節
! * 蒸気<-->雲の変換を行う.
!
! Moist adjustment.
! * Conversion from vapor to cloud.
!
xyzf_QMixAll = max( 1.0d-60, xyzf_QMixAl + xyzf_QMixBZ )
xyz_PressDry = xyz_PressDry_xyzf_xyz( xyzf_QMixAll, xyz_PressAll )
call MoistAdjustSvapPress( xyz_PressDry, xyz_ExnerNl, xyz_PTempAl, xyzf_QMixAl )
if (IdxNH4SHr /= 0) then
call MoistAdjustNH4SH( xyz_PressDry, xyz_ExnerNl, xyz_PTempAl, xyzf_QMixAl )
end if
! Set Margin
!
! call SetMargin_xyz(xyz_PTempAl)
! call SetMargin_xyzf(xyzf_QMixAl)
!-------------------------------------------
! 暖かい雨のパラメタリゼーション.
! * 蒸気<-->雨 の変換を行う
!
! Warm rain parameterization.
! * Conversion from rain to vapor.
!これまでの値を作業配列に保管
! Previous values are stored to work area.
!
xyz_PTempWork = xyz_PTempAl
xyzf_QMixWork = xyzf_QMixAl
! 雨から蒸気への混合比変化を求める
! * 温位の計算において, 混合比変化が必要となるため,
! 混合比変化を 1 つの配列として用意する.
!
! Conversion values are calculated.
!
!温度, 圧力, 混合比の全量を求める
!擾乱成分と平均成分の足し算
!
xyz_TempAll = ( xyz_PTempAl + xyz_PTempBZ ) * ( xyz_ExnerNl + xyz_ExnerBZ )
xyzf_QMixAll = max( 1.0d-60, xyzf_QMixAl + xyzf_QMixBZ )
xyz_PressDry = xyz_PressDry_xyzf_xyz( xyzf_QMixAll, xyz_PressAll )
xyzf_Rain2Gas = 0.0d0
xyzf_DelPTemp = 0.0d0
xyzf_Rain2GasNH4SH = 0.0d0
xyz_DelPTempNH4SH = 0.0d0
do s = 1, CondNum
! 値を保管
!
iG = IdxCG(s)
iC = IdxCC(s)
iR = IdxCR(s)
!飽和蒸気圧と混合比の差(飽和度)を計算.
! 雨から蒸気への変換量は飽和度に比例する.
!
xyz_NonSaturate = max( 1.0d-60, xyz_SvapPress(SpcWetID(iC), xyz_TempAll) * MolWtWet(iG) / ( MolWtDry * xyz_PressDry) - xyzf_QMixAll(:,:,:,iG) )
!雨の変換量
! 元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
!
xyzf_Rain2Gas(:,:,:,iR) = - min( DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate * ( xyzf_QMixAll(:,:,:,iR) * xyz_DensBZ )** 0.65d0, xyzf_QMixAll(:,:,:,iR) )
!蒸気の変換量
! 雨粒の変換量とは符号が逆となる
!
xyzf_Rain2Gas(:,:,:,iG) = - xyzf_Rain2Gas(:,:,:,iR)
! xyzf_DelQMix を元に潜熱を計算
!
xyzf_DelPTemp(:,:,:,s) = xyz_LatentHeat( SpcWetID(iR), xyz_TempAll ) * xyzf_Rain2Gas(:,:,:,iR) / (xyz_ExnerAll * CpDry)
end do
!飽和蒸気圧と混合比の差(飽和度)を計算.
! 雨から蒸気への変換量は飽和度に比例する.
! 未飽和度を求めたいので, マイナスをかけ算している
! (DelQMixNH4SH は, NH4SH が増加する方向, すなわち飽和度を正としている)
!
if (IdxNH4SHr /= 0) then
xyz_NonSaturate = max( 1.0d-60, - xyz_DelQMixNH4SH( xyz_TempAll, xyz_PressAll, xyz_PressDry, xyzf_QMixAll(:,:,:,IdxNH3), xyzf_QMixAll(:,:,:,IdxH2S), MolWtWet(IdxNH3), MolWtWet(IdxH2S) ) )
!雨の変換量
! 元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
!
xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) = - min( DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate * (xyzf_QMixAll(:,:,:,IdxNH4SHr) * xyz_DensBZ) ** 0.65d0, xyzf_QMixAll(:,:,:,IdxNH4SHr) )
!蒸気の変換量
! 雨粒の変換量とは符号が逆となる
!
xyzf_Rain2GasNH4SH(:,:,:,IdxNH3) = - xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) * MolWtWet(IdxNH3) / MolWtWet(IdxNH4SHr)
xyzf_Rain2GasNH4SH(:,:,:,IdxH2S) = - xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) * MolWtWet(IdxH2S) / MolWtWet(IdxNH4SHr)
xyz_DelPTempNH4SH = ReactHeatNH4SH * xyzf_Rain2GasNH4SH(:,:,:,IdxNH4SHr) / (xyz_ExnerAll * CpDry)
end if
!変化量を足し算
!
xyzf_DelQMix = xyzf_Rain2Gas + xyzf_Rain2GasNH4SH
xyz_DelPTemp = sum(xyzf_DelPTemp, 4) + xyz_DelPTempNH4SH
! 温位と混合比の計算. 雨から蒸気への変換分を追加
!
xyz_PTempAl = xyz_PTempWork + xyz_DelPTemp
xyzf_QMixAl = xyzf_QMixWork + xyzf_DelQMix
!------------------------------------------
! Output
!
xyz_PTempCond = (xyz_PTempAl - xyz_PTempOrig) / DelTime
xyzf_QMixCond = (xyzf_QMixAl - xyzf_QMixOrig) / DelTime
if ( .not. FlagDExnerDtCloud ) then
xyz_ExnerCond = 0.0d0
else
xyz_ExnerCond = xyz_DExnerDt_xyz_xyzf( xyz_PTempCond, xyzf_QMixCond )
end if
xyz_DExnerDt = xyz_DExnerDt + xyz_ExnerCond
call HistoryAutoPut(TimeN, 'PTempCond', xyz_PTempCond(1:nx, 1:ny, 1:nz))
call HistoryAutoPut(TimeN, 'ExnerCond', xyz_ExnerCond(1:nx,1:ny,1:nz))
do l = 1, ncmax
call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Cond', xyzf_QMixCond(1:nx, 1:ny, 1:nz, l))
end do
! Set Margin
!
! call SetMargin_xyz(xyz_PTempAl)
! call SetMargin_xyzf(xyzf_QMixAl)
end subroutine Cloudphys_K1969_forcing