!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005, 2006. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine WarmRainPrm
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める.
!   * 中島健介 (1994) で利用した定式をそのまま利用. 
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

module WarmRainPrm
  !
  !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める.
  !   * 中島健介 (1994) で利用した定式をそのまま利用. 
  ! 
  
  !モジュール読み込み
  use gridset, only:  DimXMin,           &! x 方向の配列の下限
    &                 DimXMax,           &! x 方向の配列の上限
    &                 DimZMin,           &! z 方向の配列の下限
    &                 DimZMax,           &! z 方向の配列の上限
    &                 SpcNum              !化学種の数
  use basicset, only: CpDry,             &!乾燥成分の比熱
    &                 MolWtWet,          &!
    &                 SpcWetID,          &!
    &                 SpcWetSymbol,      &!
    &                 xz_DensBasicZ,     &!基本場の密度
    &                 xz_PotTempBasicZ,  &!基本場の温位
    &                 xz_ExnerBasicZ,    &!基本場の無次元圧力
    &                 xz_MixRtBasicZ      !基本場の混合比
  use average,  only: xz_avr_xr        
  use differentiate_center4, only: xr_dz_xz
  use ChemData, only: ChemData_OneSpcID
  use ChemCalc, only: SvapPress, LatentHeatPerMass
  use MoistFunc,only: Vap2MixRt, xz_DelMixRtNH4SH
  
  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数を public にする
  public WarmRainPrm_Init
  public xz_LatentHeatRain2Gas
  public xz_Rain2Gas
  public xz_Cloud2Rain
  public xz_Cloud2RainNH4SH
  public xz_FallRain

  real(8)       :: MinError = 1.0d-12  !許容相対誤差
  integer      :: LoopNum      = 0
  integer      :: GasNum(10)   = 0
  integer      :: CloudNum(10) = 0
  integer      :: RainNum(10) = 0
  integer      :: NH3Num   = 0
  integer      :: H2SNum   = 0
  integer      :: NH4SHCloudNum = 0
  integer      :: NH4SHRainNum = 0
  real(8), allocatable :: RainSW(:)
  
  save MinError
  save LoopNum, RainNum, CloudNum, GasNum
  save NH3Num, H2SNum, NH4SHCloudNum, NH4SHRainNum
  save RainSW

contains  

!!!=================================================================================!!!
  subroutine WarmRainPrm_Init()

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer                  :: s
    integer                  :: n1, n2, LoopNum2

    !-----------------------------------------------------------
    ! 雨粒と雲粒と気体の ID の組を作る
    !-----------------------------------------------------------
    !初期化
    allocate( RainSW(SpcNum) )
    LoopNum  = 0
    LoopNum2 = 0
    RainSW   = 0.0d0

    !化学種の中から雨粒を作るものを選び, その配列添え字と分子量を保管.
    SelectCloud: do s = 1, SpcNum

      !'Cloud' という文字列が含まれるものの個数を数える
      n1 = index(SpcWetSymbol(s), '-Cloud' )
      if (n1 /= 0) then
        LoopNum           = LoopNum + 1
        GasNum(LoopNum)   = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID(SpcWetSymbol(s)(1:n1-3) // '-g'))
        CloudNum(LoopNum) = s
      end if

      !'Rain' という文字列が含まれるものの個数を数える
      n2 = index(SpcWetSymbol(s), '-Rain' )
      if (n2 /= 0) then
        LoopNum2          = LoopNum2 + 1
        RainNum(LoopNum2) = s
        RainSW(s)         = 1.0d0
      end if

      ! NH4SH が存在する場合は LoopNum を 1 つ減らす
      if ( trim(SpcWetSymbol(s)) == 'NH4SH-s-Cloud' ) then 
        LoopNum = LoopNum - 1
      end if

    end do SelectCloud
    
    !-----------------------------------------------------------
    ! 硫化アンモニウム, およびアンモニアと硫化水素の ID を取得
    !   'Cloud' の方が 'Rain' よりも先に存在すると仮定している
    !-----------------------------------------------------------
    NH3Num   = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH3-g'))
    H2SNum   = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('H2S-g'))
    NH4SHCloudNum = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH4SH-s'))
    NH4SHRainNum  = NH4SHCloudNum + 1

    !-----------------------------------------------------------
    ! 確認
    !-----------------------------------------------------------
    if ( LoopNum == 0 ) then 
      write(*,*) "WarmRainPrm: CloudNum = 0, please comment out of WarmRainPrm"
!      stop
    end if

    write(*,*) "WarmRainPrm_Init, LoopNum:  ", LoopNum
    write(*,*) "WarmRainPrm_Init, GasNum:   ", GasNum    
    write(*,*) "WarmRainPrm_Init, CloudNum: ", CloudNum
    write(*,*) "WarmRainPrm_Init, RainNum:  ", RainNum    
    write(*,*) "WarmRainPrm_Init, RainSW:   ", RainSW
    write(*,*) "WarmRainPrm_Init, NH3Num:   ", NH3Num
    write(*,*) "WarmRainPrm_Init, H2SNum:   ", H2SNum
    write(*,*) "WarmRainPrm_Init, NH4SHNum: ", NH4SHCloudNum
    write(*,*) "WarmRainPrm_Init, NH4SHNum: ", NH4SHRainNum

  end subroutine WarmRainPrm_Init


!!!=================================================================================!!!  
  function xz_Rain2Gas(xz_Exner, xz_PotTemp, xz_MixRt)
    !
    ! 雨粒から蒸気への変換量を計算するためのルーチン
    !
    ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で
    ! 値が正になることを保証している. また, 元々存在する以上の雨粒が
    ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている.
    !
        
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温位の擾乱成分
    real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温度の擾乱成分 
    real(8), intent(in) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分
    real(8)             :: xz_Rain2Gas(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !
    real(8)             :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温度の擾乱成分 + 平均成分
    real(8)             :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分 + 平均成分
    real(8)             :: NonSaturate    !未飽和度(飽和混合比と蒸気の混合比の差)
    integer             :: s, i, k

    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xz_TempAll  = ( xz_PotTemp + xz_PotTempBasicZ ) *  ( xz_Exner + xz_ExnerBasicZ )
    xz_MixRtAll(:,:,:)   = xz_MixRt(:,:,:) + xz_MixRtBasicZ(:,:,:)
    xz_Rain2Gas(:,:,:) = 0.0d0
    
    do s = 1, LoopNum
      do k = DimZMin, DimZmax
        do i = DimXMin, DimXMax
          
          !飽和蒸気圧と混合比の差(未飽和度)を計算. 
          !  雨から蒸気への変換量は未飽和度に応じる.
          NonSaturate =                                              &
            & Vap2MixRt(                                             &
            &     MolWtWet(CloudNum(s)),                             &
            &     SvapPress(SpcWetID(CloudNum(s)), xz_TempAll(i,k)), &
            &     (xz_Exner(i,k) + xz_ExnerBasicZ(i,k))              &
            &   )                                                    &
            & - xz_MixRtAll(i,k, GasNum(s))
  
          !雨の変換量
          !  元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
          xz_Rain2Gas(i,k,RainNum(s)) =                              &
            & - min(                                                   &
            &        4.85d-2 * max(0.0d0, NonSaturate )                &
            &        * (                                               &
            &             max( 0.0d0, xz_MixRtAll(i,k,RainNum(s)) )    &
            &             * xz_DensBasicZ(i,k)                         &
            &           ) ** 0.65d0 ,                                  &
            &        max( 0.0d0, xz_MixRtAll(i,k,RainNum(s)) )         &
            &      )

          !蒸気の変換量
          !  雨粒の変換量とは符号が逆となる
          xz_Rain2Gas(i,k,GasNum(s)) = - xz_Rain2Gas(i,k,RainNum(s)) 
        end do
      end do
    end do

!    write(*,*) 'R2V: ', minval(xz_Rain2Gas(:,:,1)), maxval(xz_Rain2Gas(:,:,1))
!    write(*,*) 'R2V: ', minval(xz_Rain2Gas(:,:,2)), maxval(xz_Rain2Gas(:,:,2))
!    write(*,*) 'R2V: ', minval(xz_Rain2Gas(:,:,3)), maxval(xz_Rain2Gas(:,:,3))
    
  end function xz_Rain2Gas
    

!!!=================================================================================!!!  
  function xz_LatentHeatRain2Gas(xz_Exner, xz_PotTemp, xz_MixRt)
    !
    ! 雨粒から蒸気への変換量を計算するためのルーチン
    !
    ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で
    ! 値が正になることを保証している. また, 元々存在する以上の雨粒が
    ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている.
    !
        
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温位の擾乱成分
    real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温度の擾乱成分 
    real(8), intent(in) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分
    real(8)             :: xz_LatentHeatRain2Gas(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !
    real(8)             :: xz_Rain2Gas(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !
    real(8)             :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !温度の擾乱成分 + 平均成分
    real(8)             :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分 + 平均成分
    real(8)             :: NonSaturate    !未飽和度(飽和混合比と蒸気の混合比の差)
    integer             :: s, i, k

    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xz_TempAll  = ( xz_PotTemp + xz_PotTempBasicZ ) *  ( xz_Exner + xz_ExnerBasicZ )
    xz_MixRtAll(:,:,:)   = xz_MixRt(:,:,:) + xz_MixRtBasicZ(:,:,:)
    xz_Rain2Gas(:,:,:) = 0.0d0
    xz_LatentHeatRain2Gas(:,:) = 0.0d0
    
    do s = 1, LoopNum
      do k = DimZMin, DimZmax
        do i = DimXMin, DimXMax
          
          !飽和蒸気圧と混合比の差(未飽和度)を計算. 
          !  雨から蒸気への変換量は未飽和度に応じる.
          NonSaturate =                                              &
            & Vap2MixRt(                                             &
            &     MolWtWet(CloudNum(s)),                             &
            &     SvapPress(SpcWetID(CloudNum(s)), xz_TempAll(i,k)), &
            &     (xz_Exner(i,k) + xz_ExnerBasicZ(i,k))              &
            &   )                                                    &
            & - xz_MixRtAll(i,k, GasNum(s))
  
          !雨の変換量
          !  元々の雨粒の混合比以上に蒸発が生じないように上限値を設定
          xz_Rain2Gas(i,k,RainNum(s)) =                              &
            & - min(                                                   &
            &        4.85d-2 * max(0.0d0, NonSaturate )                &
            &        * (                                               &
            &             max( 0.0d0, xz_MixRtAll(i,k,RainNum(s)) )    &
            &             * xz_DensBasicZ(i,k)                         &
            &           ) ** 0.65d0 ,                                  &
            &        max( 0.0d0, xz_MixRtAll(i,k,RainNum(s)) )         &
            &      )
          
          !雨から蒸気への相変化に伴う発熱
          xz_LatentHeatRain2Gas(i,k) =                                   &
            & xz_LatentHeatRain2Gas(i,k)                                 &
            & + LatentHeatPerMass( SpcWetID(RainNum(s)), xz_TempAll(i,k) ) &
            &   * xz_Rain2Gas(i,k,RainNum(s)) 
          
        end do
      end do
    end do
    
  end function xz_LatentHeatRain2Gas
  

!!!=================================================================================!!!  
  function xz_Cloud2Rain( xz_MixRt )
    !
    ! 雲粒から雨粒への変換量を計算するためのルーチン
    ! 併合成長は Kessler (1969) のパラメタリゼーションを利用し, 
    ! 衝突合体成長は Kessler (1969) のパラメタリゼーションを利用する. 
    !
    ! 変換量および, 雲粒と雨粒の混合比は正の量なので, 計算の途中途中で
    ! 値が正になることを保証している. また, 元々存在する以上の雲粒が
    ! 雨粒に変換されないように, 元々の雲粒混合比を変換量の上限としている.
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分
    real(8)             :: xz_Cloud2Rain(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !雲から雨への変換量
    real(8)             :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分 + 平均成分
    real(8)             :: xz_AutoConv(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !飽和混合比
    real(8)             :: xz_Collect(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !規格化された潜熱
!    real(8), parameter  :: N0 = 5.0d7 
!    real(8), parameter  :: D0 = 3.66d-1
    integer             :: i, k, s

    xz_MixRtAll = xz_MixRt + xz_MixRtBasicZ 
    xz_Cloud2Rain  = 0.0d0

    do s = 1, LoopNum
      xz_AutoConv = 0.0d0
      xz_Collect  = 0.0d0
      
      do k = DimZMin, DimZmax
        do i = DimXMin, DimXMax
          
          !併合成長
          !  Kessler (1969) のパラメタリゼーション        
          xz_AutoConv(i,k) =                                                &
            & 1.0d-3 * max( 0.0d0, ( xz_MixRtAll(i,k,CloudNum(s)) - 1.0d-3) )

!          !  Berry (1968) のパラメタリゼーション        
!          xz_AutoConv(i,k) =                                                &
!            & xz_DensBasicZ(i,k)                                            &
!            & * (                                                           &
!            &      max( 0.0d0, xz_MixRtAll(i,k,CloudNum(s)) ) ** 3.0d0      &
!            &    ) * 1.0d6                                                  &
!            & / 60.0d0 * (                                                  &
!            &       2.0d0 * max( 0.0d0, xz_MixRtAll(i,k,CloudNum(s)) )      &
!            &     + 60.0d0 * 2.66d-8 * N0 / ( xz_DensBasicZ(i,k) * D0 )     &
!            &    ) 
          
          !衝突合体成長
          !  Kessler (1969) のパラメタリゼーション    
          xz_Collect(i,k) =                                                        &
            &  2.2d0 * max( 0.0d0, xz_MixRtAll(i,k, CloudNum(s)) )                 &
            &  * (                                                                 &
            &       max( 0.0d0, xz_MixRtAll(i,k,RainNum(s)) ) * xz_DensBasicZ(i,k) &
            &     ) ** 0.875d0  
         
          !雲の変換量: 併合成長と合体衝突の和
          !  元々の変化量を上限値として設定する. 負の値となる.
          xz_Cloud2Rain(i,k,CloudNum(s)) =                         &
            & - min(                                               &
            &         max( 0.0d0, xz_MixRtAll(i,k,CloudNum(s)) ),  &
            &         ( xz_AutoConv(i,k) + xz_Collect(i,k) )       &
            &       )
          
          !雨の変換量. 符号は雲の変換量とは反対. 
          xz_Cloud2Rain(i,k,RainNum(s)) = - xz_Cloud2Rain(i,k,CloudNum(s)) 
          
        end do
      end do
    end do

!    write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,1)), maxval(xz_Cloud2Rain(:,:,1))
!    write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,2)), maxval(xz_Cloud2Rain(:,:,2))
!    write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,3)), maxval(xz_Cloud2Rain(:,:,3))
    
  end function xz_Cloud2Rain


!!!=================================================================================!!!
  function xz_Cloud2RainNH4SH( xz_MixRt )
    !
    ! 雲粒から雨粒への変換量を計算するためのルーチン
    ! 併合成長は Kessler (1969) のパラメタリゼーションを利用し, 
    ! 衝突合体成長は Kessler (1969) のパラメタリゼーションを利用する. 
    !
    ! 変換量および, 雲粒と雨粒の混合比は正の量なので, 計算の途中途中で
    ! 値が正になることを保証している. また, 元々存在する以上の雲粒が
    ! 雨粒に変換されないように, 元々の雲粒混合比を変換量の上限としている.
    !
    
    !暗黙の型宣言禁止
    implicit none

    !変数定義
    real(8), intent(in) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分
    real(8)             :: xz_Cloud2RainNH4SH(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !雲から雨への変換量
    real(8)             :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                          !混合比の擾乱成分 + 平均成分
    real(8)             :: xz_AutoConv(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !飽和混合比
    real(8)             :: xz_Collect(DimXMin:DimXMax, DimZMin:DimZMax)
                                          !規格化された潜熱
    integer             :: i, k

    xz_MixRtAll = xz_MixRt + xz_MixRtBasicZ 
    xz_Cloud2RainNH4SH  = 0.0d0
    xz_AutoConv = 0.0d0
    xz_Collect  = 0.0d0
      
    do k = DimZMin, DimZmax
      do i = DimXMin, DimXMax
          
        !併合成長
        !  Kessler (1969) のパラメタリゼーション        
        xz_AutoConv(i,k) =                                                &
          & 1.0d-3 * max( 0.0d0, ( xz_MixRtAll(i,k,NH4SHCloudNum) - 1.0d-3) )

        !衝突合体成長
        !  Kessler (1969) のパラメタリゼーション    
        xz_Collect(i,k) =                                                        &
          &  2.2d0 * max( 0.0d0, xz_MixRtAll(i,k, NH4SHCloudNum) )                 &
          &  * (                                                                 &
          &       max( 0.0d0, xz_MixRtAll(i,k,NH4SHRainNum) ) * xz_DensBasicZ(i,k) &
          &     ) ** 0.875d0  
        
        !雲の変換量: 併合成長と合体衝突の和
        !  元々の変化量を上限値として設定する. 負の値となる.
        xz_Cloud2RainNH4SH(i,k,NH4SHCloudNum) =                  &
          & - min(                                               &
          &         max( 0.0d0, xz_MixRtAll(i,k,NH4SHCloudNum) ),&
          &         ( xz_AutoConv(i,k) + xz_Collect(i,k) )       &
          &       )
        
        !雨の変換量. 符号は雲の変換量とは反対. 
        xz_Cloud2RainNH4SH(i,k,NH4SHRainNum) = - xz_Cloud2RainNH4SH(i,k,NH4SHCloudNum) 
        
      end do
    end do
    
!    write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,1)), maxval(xz_Cloud2Rain(:,:,1))
!    write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,2)), maxval(xz_Cloud2Rain(:,:,2))
!    write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,3)), maxval(xz_Cloud2Rain(:,:,3))
    
  end function xz_Cloud2RainNH4SH


!!!=================================================================================!!!
  function xz_FallRain( xz_MixRt, num )
    !
    ! 雨粒の落下による移流を求める. 
    ! 
    ! このルーチンの引数 xz_MixRt は 2 次元配列である. 
    ! 引数に与えられた混合比に対し, 移流を計算する. 
    ! もしも, xz_MixRt が雨粒の混合比を表すなら(RainSW = 1)なら
    ! その値を出力する. 
    ! もしも, xz_MixRt が蒸気もしくは雲粒の混合比を表すなら(RainSW = 0)
    ! ならば, 落下による移流はゼロとする. 
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    integer, intent(in) :: num
    real(8), intent(in) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax)  
                                                 !蒸気混合比(擾乱)
    real(8)             :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax)  
                                                 !蒸気混合比(擾乱 + 平均場)
    real(8)             :: xz_FallRain(DimXMin:DimXMax, DimZMin:DimZMax)
                                                 !雨粒の落下効果
    real(8)             :: xz_VelZRain(DimXMin:DimXMax, DimZMin:DimZMax)
                                                 !雨粒落下速度

    xz_MixRtAll = xz_MixRt + xz_MixRtBasicZ(:,:,num)
    xz_FallRain = 0.0d0
    xz_VelZRain = 0.0d0

    where (xz_MixRt > 1.0d-16) 
      
      !雨粒終端速度
      xz_VelZRain = 12.2d0 * (xz_MixRtAll ** 0.125d0) 
      
      !落下による移流
      !  Dens の avr を取ってから割ると, ゼロ割が生じるので注意
      xz_FallRain =                                                     &
        & - xz_avr_xr(                                                  &
        &         xr_dz_xz(xz_DensBasicZ * xz_VelZRain * xz_MixRtAll)   &
        &      ) / xz_DensBasicZ                                        &
        &        * RainSW(num)
      
    end where

!    write(*,*) 'MixRt: ', minval( xz_MixRt    ),  maxval( xz_MixRt    )
!    write(*,*) 'Fall:  ', minval( xz_FallRain ),  maxval( xz_FallRain )

  end function xz_FallRain
  
end module WarmRainPrm
