!= Module MoistBuoyancy
!
! Authors::   SUGIYAMA Ko-ichiro
! Version::   $Id: moistbuoyancy.f90,v 1.15 2006/11/22 15:40:55 sugiyama Exp $
! Tag Name::  $Name: arare4-20061224 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!
!== Error Handling
!
!== Bugs
!
!== Note
!
!
!== Future Plans
!
!

module MoistBuoyancy
  !
  !
  !

  !モジュール読み込み 
  use dc_message, only: MessageNotify

  use gridset, only:  SpcNum,            &! 化学種の数
    &                 DimXMin,           &! x 方向の配列の下限
    &                 DimXMax,           &! x 方向の配列の上限
    &                 DimZMin,           &! z 方向の配列の下限
    &                 DimZMax,           &! z 方向の配列の上限
    &                 DelX,              &! x 方向の格子点間隔
    &                 DelZ                ! z 方向の格子点間隔
  use basicset, only: MolWtDry,          &!乾燥成分の分子量
    &                 Grav,              &!重力加速度
    &                 SpcWetID,          &!凝縮成分のID
    &                 SpcWetSymbol,      &!凝縮成分の名前
    &                 MolWtWet,          &!凝縮成分の分子量
    &                 CpDry,             &!乾燥成分の平均定圧比熱 [J/K kg]
    &                 xz_PotTempBasicZ,  &!基本場の温位
    &                 xz_EffMolWtBasicZ, &!基本場の分子量効果
    &                 xz_ExnerBasicZ,    &!基本場の無次元圧力
    &                 xza_MixRtBasicZ     !基本場の混合比
  use ChemCalc, only: xz_LatentHeat       !潜熱
  use ChemData, only: ChemData_OneSpcID
  use average,  only: xr_avr_xz, xz_avr_xr
  use differentiate_center2, only: xr_dz_xz

  !暗黙の型宣言禁止
  implicit none
  
  !属性の指定
  private
  
  !関数を public に設定
  public MoistBuoy_Init
  public xz_BuoyMoistKm
  public xr_BuoyMolWt
  public xr_BuoyDrag

  !変数の定義
  integer     :: LoopNum      = 0
  integer     :: GasCount     = 0
  integer     :: GasNum(10)   = 0
  integer     :: CloudNum(10) = 0
  real(8)     :: Cm           = 2.0d-1
  real(8)     :: MixLen       = 0.0d0
  real(8), allocatable :: xz_MixRtBasicZPerMolWt(:,:)
                              !基本場の混合比 / 分子量 の総和
  real(8), allocatable :: xr_MixRtBasicZPerMolWt(:,:)
                              !基本場の混合比 / 分子量 の総和
  real(8), allocatable :: xz_MixRtBasicZ(:,:)
                              !基本場の混合比の総和
  real(8), allocatable :: xr_MixRtBasicZ(:,:)
                              !基本場の混合比の総和
  
  save GasNum, CloudNum, LoopNum
  save Cm, MixLen
  save xz_MixRtBasicZPerMolWt
  save xr_MixRtBasicZPerMolWt
  save xz_MixRtBasicZ
  save xr_MixRtBasicZ 

contains

!!!------------------------------------------------------------------!!!
  subroutine MoistBuoy_Init( )

    use ChemData, only: ChemData_OneSpcID

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer              :: s
    integer              :: n1, n2
    real(8), allocatable :: xza_MixRtBasicZPerMolWt(:,:,:)
                                  !基本場の混合比 / 分子量

    !-----------------------------------------------------------
    ! 混合距離
    !-----------------------------------------------------------
    MixLen = sqrt(DelX * DelZ) 
    
    !-----------------------------------------------------------
    ! 雲粒と気体の ID の組を作る
    !-----------------------------------------------------------
    !化学種の中から雲粒を作るものを選び, その配列添え字と分子量を保管.
    LoopNum  = 0
    GasCount = 0

    SelectCloud: do s = 1, SpcNum

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

    !-----------------------------------------------------------
    ! 確認
    !-----------------------------------------------------------
    call MessageNotify( "M", &
      & "MoistBuoy_Init", "LoopNum = %d",  i=(/LoopNum/) )
    call MessageNotify( "M", &
      & "MoistBuoy_Init", "GasCount = %d", i=(/GasCount/))
    call MessageNotify( "M", &
      & "MoistBuoy_Init", "GasNum = %d",   i=(/GasNum/)  )
    call MessageNotify( "M", &
      & "MoistBuoy_Init", "CloudNum= % ad", i=(/CloudNum/))  

    if ( LoopNum /= GasCount ) then 
      call MessageNotify( "E", &
        &                 "MoistBuoy_Init", "GasCount /= LoopNum" )
      stop
    end if

    !-----------------------------------------------------------
    ! 作業配列の初期化. 
    !-----------------------------------------------------------
    allocate( &
      & xza_MixRtBasicZPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum), &
      & xr_MixRtBasicZPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xz_MixRtBasicZPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_MixRtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax),                   &
      & xz_MixRtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)                    &
      & )

    xza_MixRtBasicZPerMolWt = 0.0d0
    xr_MixRtBasicZPerMolWt = 0.0d0
    xz_MixRtBasicZPerMolWt = 0.0d0
    xr_MixRtBasicZ = 0.0d0
    xz_MixRtBasicZ = 0.0d0

    call MessageNotify( "M", &
      & "MoistBuoy_Init", "MolWtWet = %f", d=(/pack(MolWtWet(:), .true.) /) )

    do s = 1, LoopNum
      xza_MixRtBasicZPerMolWt(:,:,GasNum(s)) = &
        & xza_MixRtBasicZ(:,:,GasNum(s)) / MolWtWet(GasNum(s))
    end do
    xz_MixRtBasicZPerMolWt = sum(xza_MixRtBasicZPerMolWt, 3) 
    xr_MixRtBasicZPerMolWt = xr_avr_xz( sum(xza_MixRtBasicZPerMolWt, 3) )
    xz_MixRtBasicZ         =  sum(xza_MixRtBasicZ, 3) 
    xr_MixRtBasicZ         = xr_avr_xz( sum(xza_MixRtBasicZ, 3) )     
    
  end subroutine MoistBuoy_Init


!!!------------------------------------------------------------------------!!!
  function xz_BuoyMoistKm(xz_PotTemp, xz_Exner, xza_MixRt)
    !
    !本来は NH4SH の凝結熱を考慮しないとならないのだが, ここでは
    !とりあえず無視する. 
    !

    !暗黙の型宣言禁止
    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)  :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                               !凝縮成分の混合比
    real(8)              :: xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                               !凝縮成分の混合比
    real(8)              :: xz_BuoyMoistKm(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !
    real(8)              :: xza_LatentHeat(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum)
                                               !潜熱
    real(8)              :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !温度
    real(8)              :: xz_EffHeat(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !
    real(8)              :: xz_EffPotTemp(DimXMin:DimXMax, DimZMin:DimZMax)    
    real(8)              :: xz_EffMolWt(DimXMin:DimXMax, DimZMin:DimZMax)    
                                               !
    real(8)              :: xza_MixRtPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum)
                                               !混合比/分子量
    integer              :: s
    
    
    !温度, 圧力, 混合比の全量を求める
    !擾乱成分と平均成分の足し算
    xz_TempAll = ( xz_PotTemp + xz_PotTempBasicZ ) * ( xz_Exner + xz_ExnerBasicZ )
    xza_MixRtAll = xza_MixRtBasicZ + xza_MixRt
    xza_LatentHeat = 0.0d0
    
    !作業配列の初期化. 気体のみ利用
    do s = 1, LoopNum
      xza_MixRtPerMolWt(:,:,GasNum(s)) =                 &
        & xza_MixRt(:,:,GasNum(s)) / MolWtWet(GasNum(s))
    end do
    
    !温度の効果
    xz_EffPotTemp = xz_PotTemp / xz_PotTempBasicZ 
    
    !分子量効果 + 引きづりの効果
    xz_EffMolWt =                                           &
      & + sum(xza_MixRtPerMolWt, 3)                         &
      &    / ( 1.0d0 / MolWtDry + xz_MixRtBasicZPerMolWt )  &
      & - sum(xza_MixRt, 3) / ( 1.0d0 + xz_MixRtBasicZ )
    
    !蒸気が蒸発する場合の潜熱を計算
    !  分子量の部分はいつでも効くが潜熱は飽和していないと効かないので, 
    !  雲の混合比がゼロの時には, 潜熱の寄与はゼロとなるように調節している
    do s = 1, LoopNum
      xza_LatentHeat(:,:,s) =                                 &
        & xz_LatentHeat( SpcWetID(CloudNum(s)), xz_TempAll )  &
        &  * xza_MixRtAll(:,:,GasNum(s))                      &
        &  * ( 5.0d-1 + sign( 5.0d-1, (xza_MixRtAll(:,:,CloudNum(s)) - 1.0d-4) ) )
    end do
    xz_EffHeat = sum( xza_LatentHeat, 3 ) * xz_EffMolWtBasicZ &
      &           / ( CpDry * xz_ExnerBasicZ ) 
    
    !乱流拡散係数の時間発展式の浮力項を決める
    xz_BuoyMoistKm = &
      &  - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 )     &
      &    * xz_avr_xr(                                             &
      &        xr_dz_xz(                                            &
      &           xz_EffHeat                                        &
      &         + xz_PotTempBasicZ / xz_EffMolWtBasicZ              &
      &           * ( 1.0d0 + xz_EffPotTemp + xz_EffMolWt )         &
      &         )                                                   &
      &      )                                                      &
      &    / ( 2.0d0 * xz_PotTempBasicZ / xz_EffMolWtBasicZ)                   

  end function xz_BuoyMoistKm
    

!!!------------------------------------------------------------------------!!!
  function xr_BuoyMolWt(xza_MixRt)
    !
    ! 鉛直方向の運動方程式に現れる浮力項のうち, 
    ! 分子量の効果だけを求める
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in)  :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                               !凝縮成分の混合比
    real(8)              :: xr_BuoyMolWt(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !浮力項(分子量効果)
    real(8)              :: xza_MixRtPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum)
                                               !混合比/分子量
    integer              :: s
    
    !初期化
    xr_BuoyMolWt           = 0.0d0

    !作業配列の初期化. 気体のみ利用
    do s = 1, LoopNum
      xza_MixRtPerMolWt(:,:,GasNum(s)) =                 &
        & xza_MixRt(:,:,GasNum(s)) / MolWtWet(GasNum(s))
    end do
    
    !浮力項の計算
    xr_BuoyMolWt =                                             &
      & + Grav * xr_avr_xz( sum(xza_MixRtPerMolWt, 3) )        &
      &    / ( 1.0d0 / MolWtDry + xr_MixRtBasicZPerMolWt )     &
      & - Grav * xr_avr_xz( sum(xza_MixRt(:,:,1:LoopNum), 3) ) &
      &    / ( 1.0d0 + xr_MixRtBasicZ )

  end function xr_BuoyMolWt


!!!------------------------------------------------------------------------!!!
  function xr_BuoyDrag(xza_MixRt) 
    !
    ! 鉛直方向の運動方程式に現れる浮力項のうち, 
    ! 引きづりのの効果だけを求める
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8), intent(in)  :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                               !凝縮成分の混合比
    real(8)              :: xr_BuoyDrag(DimXMin:DimXMax, DimZMin:DimZMax)
                                               !浮力項(分子量効果)
    !初期化
    xr_BuoyDrag           = 0.0d0

    !浮力項の計算
    xr_BuoyDrag =                                                     &
      & - Grav * xr_avr_xz( sum(xza_MixRt(:,:,LoopNum+1:SpcNum), 3) ) &
      &    / ( 1.0d0 + xr_MixRtBasicZ )

  end function xr_BuoyDrag

end module MoistBuoyancy
