!= Module cloudphys_k1969
!
! Authors::   杉山耕一朗(SUGIYAMA Ko-ichiro), 小高正嗣 (ODAKA Masatsugu), 高橋芳幸 (YOSHIYUKI Takahashi)
! Version::   $Id: cloudphys_k1969.f90,v 1.22 2012-12-19 05:15:02 odakker Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める.
!   * 中島健介 (1994) で利用した定式をそのまま利用.
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

module cloudphys_k1969
  !
  !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める.
  !   * 中島健介 (1994) で利用した定式をそのまま利用.
  !

  !モジュール読み込み
  use dc_types,   only : DP, STRING
  use dc_iounit,  only : FileOpen
  use dc_message, only : MessageNotify
  use gtool_historyauto, only: HistoryAutoAddVariable, HistoryAutoPut

  use mpi_wrapper,only: myrank
  use timeset, only:  DelTimeLong, TimeN
  use gridset, only : imin,              &!x 方向の配列の下限
    &                 imax,              &!x 方向の配列の上限
    &                 jmin,              &!y 方向の配列の上限
    &                 jmax,              &!y 方向の配列の上限
    &                 kmin,              &!z 方向の配列の下限
    &                 kmax,              &!z 方向の配列の上限
    &                 nx, ny, nz, ncmax      !物理領域の大きさ
  use constants,only: PressBasis,        &!温位の基準圧力
    &                 CpDry,             &!乾燥成分の比熱
    &                 MolWtDry,          &!
    &                 GasRDry             !乾燥成分の気体定数
  use basicset, only: xyz_DensBZ,        &!基本場の密度
    &                 xyz_PTempBZ,       &!基本場の温位
    &                 xyz_ExnerBZ,       &!基本場の無次元圧力
    &                 xyzf_QMixBZ,       &!基本場の混合比
    &                 xyz_VPTempBZ        !温位の基本場
  use composition, only:                    &
    &                 MolWtWet,          &!
    &                 SpcWetID,          &!
    &                 SpcWetSymbol,      &!
    &                 GasNum,            &!
    &                 CondNum,           &!凝結過程の数
    &                 IdxG,              &!蒸気の配列添え字
    &                 IdxCG,             &!凝結過程(蒸気)の配列添え字
    &                 IdxCC,             &!凝結過程(雲)の配列添え字
    &                 IdxCR,             &!凝結過程(雲)の配列添え字
    &                 CloudNum,          &!雲の数
    &                 RainNum,           &!雨の数
    &                 IdxC,              &!雲の配列添え字
    &                 IdxR,              &!雨の配列添え字
    &                 IdxNH3,            &!NH3(蒸気)の配列添え字
    &                 IdxH2S,            &!H2S(蒸気)の配列添え字
    &                 IdxNH4SHr           !NH4SH(雨)の配列添え字
  use axesset, only : xyr_avr_xyz
  use xyz_deriv_module,only : xyz_dz_xyr
  use ChemCalc,  only : xyz_SvapPress, xyz_LatentHeat, ReactHeatNH4SH, xyz_DelQMixNH4SH
  use MoistAdjust, only: MoistAdjustSvapPress, MoistAdjustNH4SH
  use namelist_util, only: namelist_filename
  use DExnerDt, only: xyz_DExnerDt_xyzf, xyz_DExnerDt_xyz_xyzf

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public にする
  public Cloudphys_K1969_Init
  public Cloudphys_K1969_forcing
  public Cloudphys_K1969_FallRain

  real(DP), save :: FactorJ      = 1.0d0 !雲物理過程のパラメータ
                                         !木星では 3.0d0
                                         !地球では 1.0d0 とする
  real(DP), save :: AutoConvTime = 1.0d3 !併合成長の時定数 [sec]
  real(DP), save :: QMixCr       = 1.0d-3
                                         !併合成長を生じる臨界混合比 [kg/kg]
  logical, save :: FlagDExnerDtCloud = .true.
  logical, save :: FlagDExnerDtFall  = .true.

contains

!!!=================================================================================!!!
  subroutine Cloudphys_K1969_Init

    !暗黙の型宣言禁止
    implicit none

    !内部変数
    integer  :: unit    !装置番号
    integer  :: l
    character(STRING) :: Planet = ""
    character(*), parameter:: module_name = 'Cloudphys_K1969_Init'

    !-----------------------------------------------------------
    ! NAMELIST から情報を取得
    !-----------------------------------------------------------
    ! NAMELIST から情報を取得
    !
    NAMELIST /cloudphys_k1969_nml/ &
      & Planet, FactorJ, AutoConvTime, QMixCr, &
      & FlagDExnerDtCloud, FlagDExnerDtFall

    call FileOpen(unit, file=namelist_filename, mode='r')
    read(unit, NML=cloudphys_k1969_nml)
    close(unit)

    if (trim(Planet) == "Earth") then
      FactorJ = 1.0d0
    elseif (trim(Planet) == "Jupiter") then
      FactorJ = 3.0d0
    end if

    ! 出力
    !
    if (myrank == 0) then
      call MessageNotify( "M", &
        &  module_name, "Planet = %c",  c1=trim(Planet))
      call MessageNotify( "M", &
        &  module_name, "FactorJ = %f",  d=(/FactorJ/) )
      call MessageNotify( "M", &
        &  module_name, "AutoConvTime = %f",  d=(/AutoConvTime/) )
      call MessageNotify( "M", &
        &  module_name, "QMixCr = %f",  d=(/QMixCr/) )
      call MessageNotify( "M", &
        & module_name, "FlagDExnerDtCloud= %b", l=(/ FlagDExnerDtCloud /))
      call MessageNotify( "M", &
        & module_name, "FlagDExnerDtFall= %b", l=(/ FlagDExnerDtFall /))
     end if

    ! ヒストリデータ定義
    !
    call HistoryAutoAddVariable(  &
      & varname='ExnerFall',      &
      & dims=(/'x','y','z','t'/), &
      & longname='Fall term of Exner function', &
      & units='kg.kg-1.s-1',            &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='PTempCond',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Latent heat term of potential temperature', &
      & units='K.s-1',    &
      & xtype='float')

    call HistoryAutoAddVariable(  &
      & varname='ExnerCond',&
      & dims=(/'x','y','z','t'/),     &
      & longname='Latent heat term of exner function', &
      & units='K.s-1',    &
      & xtype='float')

    do l = 1, ncmax
      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_Cond', &
        & dims=(/'x','y','z','t'/),     &
        & longname='Condensation term of '          &
        &           //trim(SpcWetSymbol(l))//' mixing ratio',  &
        & units='kg.kg-1.s-1',    &
        & xtype='float')

      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_Fall', &
        & dims=(/'x','y','z','t'/),     &
        & longname='Fall Rain term of '          &
        &           //trim(SpcWetSymbol(l))//' mixing ratio',  &
        & units='kg.kg-1.s-1',    &
        & xtype='float')

      call HistoryAutoAddVariable(  &
        & varname=trim(SpcWetSymbol(l))//'_FallFluxAtLB', &
        & dims=(/'x','y','t'/),     &
        & longname='Falling Rain Flux '          &
        &           //trim(SpcWetSymbol(l)),  &
        & units='kg.m-2.s-1',    &
        & xtype='float')
    end do

  end subroutine Cloudphys_K1969_Init
!!!=================================================================================!!!

  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,            & ! (in)
      & xyz_ExnerNl,             & ! (in)
      & xyz_PTempAl,             & ! (inout)
      & xyzf_QMixAl              & ! (inout)
      & )
    if (IdxNH4SHr /= 0) then
      call MoistAdjustNH4SH(     &
        & xyz_PressDry,          & !(in)
        & xyz_ExnerNl,           & !(in)
        & xyz_PTempAl,           & !(inout)
        & xyzf_QMixAl            & !(inout)
        & )
    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 ( 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

    ! Set Margin
    !
!    call SetMargin_xyz(xyz_PTempAl)
!    call SetMargin_xyzf(xyzf_QMixAl)

  end subroutine Cloudphys_K1969_forcing

!!!=================================================================================!!!
  subroutine CloudPhys_K1969_FallRain(xyzf_QMix, xyz_DExnerDt, xyzf_DQMixDt )
    !
    ! 雨粒の落下による移流を求める.
    !

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(DP), intent(in) :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                                 !蒸気混合比(擾乱)
    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)  :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax, ncmax)
                                                 !蒸気混合比(擾乱 + 平均場)
    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)
    integer  :: s, l

    xyzf_QMixAll = max( 0.0d0, xyzf_QMix + xyzf_QMixBZ )

    xyrf_FallRainFlux = 0.0d0
    do s = 1, RainNum
      ! 雨粒終端速度
      xyz_VelZRain = - 12.2d0 * FactorJ               &
        * ( xyzf_QMixAll(:,:,:,IdxR(s)) ** 0.125d0 )

      xyrf_FallRainFlux(:,:,:,IdxR(s)) =              &
        &  xyr_avr_xyz(                               &
        &               xyz_DensBZ                    &
        &               * xyzf_QMixAll(:,:,:,IdxR(s)) &
        &               * xyz_VelZRain                &
        &              )
    end do
    ! 上端のフラックスはゼロ
    do s = 1, RainNum
      xyrf_FallRainFlux(:,:,nz,IdxR(s)) = 0.0d0
    end do

    ! 雨粒落下による時間変化率
    xyzf_FallRain = 0.0d0
    do s = 1, RainNum
      xyzf_FallRain(:,:,:,IdxR(s)) =                       &
        & - xyz_dz_xyr( xyrf_FallRainFlux(:,:,:,IdxR(s)) ) &
        & / xyz_DensBZ
    end do

    xyzf_DQMixDt = xyzf_DQMixDt + xyzf_FallRain

    if ( FlagDExnerDtFall ) then
      xyz_ExnerFall = xyz_DExnerDt_xyzf( xyzf_DQMixDt )
    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

    ! SetMargin
    !
!    call SetMargin_xyzf(xyzf_DQMixDt)

  end subroutine CloudPhys_K1969_FallRain

!!!----------------------------------------------------------------
  function xyz_PressDry_xyzf_xyz( xyzf_QMixAll, xyz_PressAll )
    ! 乾燥成分の分圧の計算
    !
    implicit none

    real(DP), intent(in) :: xyzf_QMixAll(imin:imax,jmin:jmax,kmin:kmax,ncmax)
    real(DP), intent(in) :: xyz_PressAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyz_PressDry_xyzf_xyz(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)             :: xyzf_QMixAllPerMolWt(imin:imax,jmin:jmax,kmin:kmax,ncmax)
    integer              :: f, iG

    ! 混合比/分子量を計算
    !
    do f = 1, GasNum
      iG = IdxG(f)
      xyzf_QMixAllPerMolWt(:,:,:,f) = xyzf_QMixAll(:,:,:,iG) / MolWtWet(iG)
    end do

    ! 乾燥成分の分圧.
    !
    xyz_PressDry_xyzf_xyz = xyz_PressAll / (1.0d0 + MolWtDry * sum( xyzf_QMixAllPerMolWt(:,:,:,1:GasNum), 4))

  end function xyz_PressDry_xyzf_xyz

end module Cloudphys_k1969