subroutine Cloudphys_K1969_forcing(xyz_ExnerNl, xyzf_QMixNl, xyz_DExnerDt, xyz_PTempAl, xyzf_QMixAl)
implicit none
real(DP), intent(in) :: xyz_ExnerNl(imin:imax, jmin:jmax, kmin:kmax)
real(DP), intent(in) :: xyzf_QMixNl(imin:imax, jmin:jmax, kmin:kmax, ncmax)
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)
real(DP) :: xyzf_FallRain(imin:imax,jmin:jmax,kmin:kmax, ncmax)
!雨粒の落下効果
real(DP) :: xyz_VelZRain(imin:imax,jmin:jmax,kmin:kmax)
!雨粒落下速度
real(DP) :: xyrf_FallRainFlux(imin:imax,jmin:jmax,kmin:kmax, ncmax)
!雨粒落下フラックス
real(DP) :: xyz_ExnerFall(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyz_QMixSat(imin:imax,jmin:jmax,kmin:kmax)
real(DP) :: xyz_QMixHum(imin:imax,jmin:jmax,kmin:kmax)
!-------------------------------------------------------------
! 時間刻み幅. Leap-frog なので, 2 \del t
!
DelTime = 2.0d0 * DelTimeLong
!-------------------------------------------------------------
! 雨粒の落下による混合比の時間変化を計算
! QMixNl を用いて QMixAl を更新する.
!
! 初期化
! 落下する分としては, Nl の値を利用する.
!
xyzf_QMixAll = max( 0.0d0, xyzf_QMixNl + xyzf_QMixBZ )
xyrf_FallRainFlux = 0.0d0
xyzf_FallRain = 0.0d0
xyzf_QMixWork = xyzf_QMixAl + xyzf_QMixBZ
do s = 1, RainNum
iR = IdxR(s)
! 雨粒終端速度
xyz_VelZRain = - 12.2d0 * FactorJ * ( xyzf_QMixAll(:,:,:,iR) ** 0.125d0 )
! フラックスの計算
!
xyrf_FallRainFlux(:,:,:,iR) = xyr_avr_xyz( xyz_DensBZ * xyz_VelZRain * xyzf_QMixAll(:,:,:,iR) )
! 上端のフラックスはゼロ
!
xyrf_FallRainFlux(:,:,nz,iR) = 0.0d0
! 雨粒落下による時間変化 (DelTime をかけてある)
!
xyzf_FallRain(:,:,:,iR) = - xyz_dz_xyr( xyrf_FallRainFlux(:,:,:,iR) ) / xyz_DensBZ * DelTime
! 元々の存在量より減ることは無いようにする.
!
xyzf_FallRain(:,:,:,iR) = max( - xyzf_QMixWork(:,:,:,iR), xyzf_FallRain(:,:,:,iR) )
end do
! 落下を考慮することで, Al の値を更新
!
xyzf_QMixAl = xyzf_QMixAl + xyzf_FallRain
! 圧力に対する tendency を求める
!
if ( FlagDExnerDtFall ) then
xyz_ExnerFall = xyz_DExnerDt_xyzf( xyzf_FallRain )
else
xyz_ExnerFall = 0.0d0
end if
! 保管
!
xyz_DExnerDt = xyz_DExnerDt + xyz_ExnerFall
call HistoryAutoPut(TimeN, 'ExnerFall', xyz_ExnerFall(1:nx,1:ny,1:nz))
do l = 1, ncmax
call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_Fall', xyzf_FallRain(1:nx, 1:ny, 1:nz, l))
call HistoryAutoPut(TimeN, trim(SpcWetSymbol(l))//'_FallFluxAtLB', xyrf_FallRainFlux(1:nx, 1:ny, 0, l))
end do
!------------------------------------------
! 初期値を保管 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( 0.0d0, xyzf_QMixAl + xyzf_QMixBZ )
xyzf_Cloud2Rain = 0.0d0
do s = 1, CloudNum
! 値を保管
!
iC = IdxC(s)
iR = IdxR(s)
!併合成長
!
xyz_AutoConv = DelTime / AutoConvTime * max( 0.0d0, ( 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)
!-------------------------------------------
! 暖かい雨のパラメタリゼーション.
! * 蒸気<-->雨 の変換を行う
!
! 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( 0.0d0, 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( 0.0d0, 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( 0.0d0, - 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
!-------------------------------------------
! 湿潤飽和調節
! * 蒸気<-->雲の変換を行う.
!
! Moist adjustment.
! * Conversion from vapor to cloud.
!
xyzf_QMixAll = max( 0.0d0, 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)
!------------------------------------------
! Output
!
xyz_PTempCond = (xyz_PTempAl - xyz_PTempOrig) / DelTime
xyzf_QMixCond = (xyzf_QMixAl - xyzf_QMixOrig) / DelTime
if ( FlagDExnerDtCloud ) then
xyz_ExnerCond = xyz_DExnerDt_xyz_xyzf( xyz_PTempCond, xyzf_QMixCond )
else
xyz_ExnerCond = 0.0d0
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
!----------------------------------------------------------------
! 飽和蒸気圧と平衡定数
!----------------------------------------------------------------
xyz_TempAll = ( xyz_PTempAl + xyz_PTempBZ ) * ( xyz_ExnerNl + xyz_ExnerBZ )
xyzf_QMixAll = max( 0.0d0, xyzf_QMixAl + xyzf_QMixBZ )
xyz_PressDry = xyz_PressDry_xyzf_xyz( xyzf_QMixAll, xyz_PressAll )
do s = 1, CondNum
iG = IdxCG(s)
iC = IdxCC(s)
!飽和蒸気圧
xyz_QMixSat = xyz_SvapPress(SpcWetID(iC), xyz_TempAll) * MolWtWet(iC) / MolWtDry / xyz_PressDry
!相対湿度
xyz_QMixHum = xyzf_QMixAll(:,:,:,iG) / xyz_QMixSat * 100.0d0
!出力
call HistoryAutoPut(TimeN, trim(SpcWetSymbol(s))//'_Hum', xyz_QMixHum(1:nx, 1:ny, 1:nz))
end do
! Set Margin
!
! call SetMargin_xyz(xyz_PTempAl)
! call SetMargin_xyzf(xyzf_QMixAl)
end subroutine Cloudphys_K1969_forcing