!= Module ECCM
!
! Authors::   SUGIYAMA Koichiro, ODAKA Masatsugu
! Version::   $Id: eccm.f90,v 1.5 2011-06-03 06:11:22 sugiyama Exp $
! Tag Name::  $Name: arare5-20111010 $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview 
!
!断熱的に上昇する気塊の温度減率を計算し, 静水圧平衡から圧力を決める
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * 比熱は定数, 平均分子量は配列
!  * オイラースキームだと精度が足りないので, ルンゲクッタスキームを利用. 
!
!== Future Plans
!

module ECCM

  !モジュール読み込み
  use dc_types, only : DP
  use dc_trace, only : DbgMessage
  use dc_message, only: MessageNotify

  use gridset,only: kmin,       &! 配列の Z 方向の下限
    &                 kmax,       &! 配列の Z 方向の上限 
    &                 ncmax
  use axesset, only:  z_dz, r_dz      !
  use constants,only: MolWtDry,      &!
    &                 CpDryMol,      &!
    &                 TempSfc,       &!
    &                 PressSfc,      &!
    &                 Grav            !
  use chemcalc, only: SvapPress,        &!
    &                 LatentHeatPerMol, &!
    &                 ReactHeatNH4SHPerMol, &
    &                 DelMolFrNH4SH 
  use composition, only: CondNum,          &!凝結過程の数
    &                 SpcWetID,         &!
    &                 IdxCG,            &!凝結過程(蒸気)の配列添え字
    &                 IdxCC,            &!凝結過程(雲)の配列添え字
    &                 GasNum,           &!気体の数
    &                 IdxNH3,           &!NH3(蒸気)の配列添え字
    &                 IdxH2S             !H2S(蒸気)の配列添え字
  use ChemData, only: GasRUniv        

  !暗黙の型宣言禁止
  implicit none

  !属性の指定
  private

  !関数の公開
  public ECCM_MolFr
  public ECCM_Stab
  public ECCM_Dry
  public ECCM_Wet

contains

!!!------------------------------------------------------------------------------!!!
  subroutine ECCM_Dry( a_MolFrIni, Humidity, z_Temp, z_Press, za_MolFr )
    !
    !== 概要
    !  * 乾燥断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める
    !  * 比熱は乾燥気塊のもので代表させる
    !    * 流体の方程式において比熱は乾燥成分のもので代表させているため
    !  * 大気の平均分子量には湿潤成分の分子量を効果
    !    * 流体の方程式において, 湿潤成分の分子量は考慮しているため
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    real(DP), intent(in) :: a_MolFrIni(1:ncmax)    !下部境界でのモル比
    real(DP), intent(in) :: Humidity                !相対湿度 ( Humidity <= 1.0 )
    real(DP), intent(out):: z_Temp(kmin:kmax) !温度
    real(DP), intent(out):: z_Press(kmin:kmax)!圧力
    real(DP), intent(out):: za_MolFr(kmin:kmax, 1:ncmax) 
                                                   !モル分率
    real(DP)             :: SatPress                !飽和蒸気圧
    real(DP)             :: VapPress                !蒸気圧
    real(DP)             :: DelMolFr
    integer             :: k, s

    !-------------------------------------------------------------
    ! 配列の初期化
    !-------------------------------------------------------------
    !地表面の分子量を決める
    za_MolFr  = 1.0d-60
    za_MolFr(1, 1:ncmax) = a_MolFrIni(1:ncmax) 

    !地表面での温度(nx は, 高度 DelZ / 2 に相当)
    z_Temp    = 1.0d-60
    z_Temp(1) = TempSfc - Grav * MolWtDry &
      &               / CpDryMol * ( z_dz(1) * 5.0d-1 )

!    write(*,*) z_Temp(1)
    
    !地表面での圧力(1 は, 高度 DelZ / 2 に相当)
    z_Press    = 1.0d-60
    z_Press(1) = PressSfc *((TempSfc / z_Temp(1)) ** (- CpDryMol /  GasRUniv))

!    write(*,*) z_Press(1)
    
    !-----------------------------------------------------------
    ! (1) 乾燥断熱線に沿った温度を決める
    ! (2) 静水圧平衡から圧力を求める
    ! (3) (1),(2) の温度圧力に対して, とある相対湿度となるモル比を決める
    !-----------------------------------------------------------    
    DtDz: do k = 1, kmax-1

      !(1)乾燥断熱線に沿って k+1 での温度を計算
      z_Temp(k+1) = z_Temp(k) - Grav * MolWtDry / CpDryMol * r_dz(k)
      
      !念為
      if (z_Temp(k+1) <= 0.0d0 ) z_Temp(k+1) = z_Temp(k) 
      
      !(2)圧力を静水圧平衡から計算
      z_Press(k+1) =                                                  &
        &  z_Press(k) * ((z_Temp(k) / z_Temp(k+1)) ** (- CpDryMol / GasRUniv)) 

      !(3)モル比の計算
      !  まずはモル比は変化しないものとしてモル比を与える
      !  飽和蒸気圧と平衡定数との平衡条件の前に適用しておく
      za_MolFr(k+1,:) = za_MolFr(k,:)
      
      do s = 1, CondNum      
        !飽和蒸気圧
        SatPress = SvapPress( SpcWetID(IdxCC(s)), z_Temp(k+1) )        
        
        !元々のモル分率を用いて現在の蒸気圧を計算
        VapPress = za_MolFr(k,IdxCG(s)) * z_Press(k+1)
        
        !飽和蒸気圧と圧力から現在のモル比を計算
        if ( VapPress > SatPress ) then         
          za_MolFr(k+1,IdxCG(s)) = max(SatPress * Humidity / z_Press(k+1), 1.0d-16)
        end if
      end do
      
      !NH4SH の平衡条件
      if ( IdxNH3 /= 0 ) then 
        DelMolFr =                                              &
          & max (                                               &
          &    DelMolFrNH4SH(                                   &
          &         z_Temp(k+1), z_Press(k+1),                  &
          &         za_MolFr(k+1,IdxNH3), za_MolFr(k+1,IdxH2S), &
          &         Humidity                                    &
          &      ),                                             &
          &    0.0d0                                            &
          &  )
        za_MolFr(k+1,IdxNH3) = za_MolFr(k+1,IdxNH3) - DelMolFr
        za_MolFr(k+1,IdxH2S) = za_MolFr(k+1,IdxH2S) - DelMolFr
      end if
      
    end do DtDz
    
  end subroutine ECCM_Dry


!!------------------------------------------------------------------------------!!!
  subroutine ECCM_Wet( a_MolFrIni, Humidity, z_Temp, z_Press, za_MolFr )
    !
    !== 概要
    !  * 湿潤断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める
    !  * 比熱は乾燥気塊のもので代表させる
    !    * 流体の方程式において比熱は乾燥成分のもので代表させているため
    !  * 大気の平均分子量には湿潤成分の分子量を効果
    !    * 流体の方程式において, 湿潤成分の分子量は考慮しているため
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    real(DP), intent(in) :: a_MolFrIni(1:ncmax)    !下部境界でのモル比
    real(DP), intent(in) :: Humidity                !相対湿度 ( Humidity <= 1.0 )
    real(DP), intent(out):: z_Temp(kmin:kmax) !温度
    real(DP), intent(out):: z_Press(kmin:kmax)!圧力
                                                   !平均分子量
    real(DP), intent(out):: za_MolFr(kmin:kmax, 1:ncmax) 
                                                   !モル分率
    real(DP)             :: SatPress                !飽和蒸気圧
    real(DP)             :: VapPress                !蒸気圧
    real(DP)             :: DelMolFr
    real(DP)             :: a_MolFr(ncmax)         !モル比の作業配列
    integer             :: k, s

    real(DP)             :: Temp1, Press1, DTempDZ1
    real(DP)             :: Temp2, Press2, DTempDZ2
    real(DP)             :: Temp3, Press3, DTempDZ3
    real(DP)             :: Temp4, Press4, DTempDZ4
    real(DP)             :: DTempDZ
    
    !-------------------------------------------------------------
    ! 配列の初期化
    !-------------------------------------------------------------
    !地表面の分子量を決める
    za_MolFr  = 1.0d-60
    za_MolFr(1, 1:ncmax)   = a_MolFrIni(1:ncmax) 
    
    !地表面での温度(1 は, 高度 DelZ / 2 に相当)
    z_Temp          = 1.0d-60
    z_Temp(1) = TempSfc - Grav * MolWtDry  &
      &               / CpDryMol * ( z_dz(1) * 5.0d-1 )
    
    !地表面での圧力(1 は, 高度 DelZ / 2 に相当)
    z_Press           = 1.0d-60
    z_Press(1)  = &
      & PressSfc *((TempSfc / z_Temp(1)) ** (- CpDryMol /  GasRUniv))
    
    !-----------------------------------------------------------
    ! 断熱減率 dT/dz の計算. 
    !-----------------------------------------------------------    
    DtDz: do k = 1, kmax-1
      
      !初期化
      za_MolFr(k+1,:) = za_MolFr(k,:)
      
      !----------------------------------------------------
      !湿潤断熱減率をルンゲクッタ法を用いて計算
      !----------------------------------------------------
      ! (0) 高度 k での値を作業配列に保管
      Temp1  = z_Temp(k)
      Press1 = z_Press(k)
      a_MolFr  = za_MolFr(k,:)

      ! (1) 高度 k での値を用いて温度変化を計算
      call ECCM_DTempDZ( Temp1, Press1, r_dz(k), a_MolFr, DTempDZ1 )

      ! (2) (1) で求めた値を用いて, 高度 k + Δk/2 での値を用いて温度変化を計算
      !     このとき, 分子量は変化しないものとする.
      Temp2  = Temp1 + DTempDZ1 * r_dz(k) * 5.0d-1
      Press2 =                                     &
        & Press1 * ((Temp1 / Temp3) ** (Grav * MolWtDry / (GasRUniv * DTempDZ2))) 
      call ECCM_DTempDZ( Temp2, Press2, r_dz(k), a_MolFr, DTempDZ2 )

      ! (3) (2) で求めた値を用いて, 高度 k + Δk/2 での値を用いて温度変化を計算
      !     このとき, 分子量は変化しないものとする.
      Temp3  = Temp1 + DTempDZ2 * r_dz(k) * 5.0d-1
      Press3 =                                                                   &
        & Press1 * ((Temp1 / Temp3) ** (Grav * MolWtDry / (GasRUniv * DTempDZ2)))
      call ECCM_DTempDZ( Temp3, Press3, r_dz(k), a_MolFr, DTempDZ3 )
      
      ! (4) (3) で求めた値を用いて, 高度 k + Δk での値を用いて温度変化を計算
      !     このとき, 分子量は変化しないものとする.
      Temp4  = Temp1 + DTempDZ3 * r_dz(k)
      Press4 =                                               &
        & Press1 * ((Temp1 / Temp4) ** (Grav * MolWtDry / (GasRUniv * DTempDZ3)))
      call ECCM_DTempDZ( Temp4, Press4, r_dz(k), a_MolFr, DTempDZ4 )
      
      ! (5) 最終的な傾きを求める
      DTempDZ = (DTempDZ1 + DTempDZ2 * 2.0d0 + DTempDZ3 * 2.0d0 + DTempDZ4) / 6.0d0

      !----------------------------------------------------
      !得られた温度減率より温度と圧力を決める
      !----------------------------------------------------
      !温度を計算
      z_Temp(k+1) = z_Temp(k) + DTempDz * r_dz(k)

      !為念
      if(z_Temp(k+1) < 0.0d0) z_Temp(k+1) = z_Temp(k) 
      
      !圧力を静水圧平衡から計算
      z_Press(k+1) =                                                  &
        &  z_Press(k) * ( ( z_Temp(k) / z_Temp(k+1))                  &
        &    ** (Grav * MolWtDry / ( DTempDZ * GasRUniv ) ) )
      
      !----------------------------------------------------
      !モル比の計算
      !----------------------------------------------------
      do s = 1, CondNum      
        !飽和蒸気圧
        SatPress = SvapPress( SpcWetID(IdxCC(s)), z_Temp(k+1) )
        
        !元々のモル分率を用いて現在の蒸気圧を計算
        VapPress = za_MolFr(k,IdxCG(s)) * z_Press(k+1)
        
        !飽和蒸気圧と圧力から現在のモル比を計算
        if ( VapPress > SatPress ) then         
          za_MolFr(k+1,IdxCG(s)) = max(SatPress * Humidity / z_Press(k+1), 1.0d-16)
        end if
      end do
      
      !NH4SH の平衡条件
      if ( IdxNH3 /= 0 ) then 
        DelMolFr =                                              &
          & max (                                               &
          &    DelMolFrNH4SH(                                   &
          &         z_Temp(k+1), z_Press(k+1),                  &
          &         za_MolFr(k+1,IdxNH3), za_MolFr(k+1,IdxH2S),       &
          &         Humidity                                    &
          &      ),                                             &
          &    0.0d0                                            &
          &  )
        za_MolFr(k+1,IdxNH3) = za_MolFr(k+1,IdxNH3) - DelMolFr
        za_MolFr(k+1,IdxH2S) = za_MolFr(k+1,IdxH2S) - DelMolFr
      end if
      
    end do DtDz
    
  end subroutine ECCM_Wet

!!!------------------------------------------------------------------------------!!!
  subroutine ECCM_MolFr( a_MolFrIni, Humidity, z_Temp, z_Press, za_MolFr )
    !
    ! 与えられた温度に対し, 気塊が断熱的に上昇した時に実現される
    ! モル比のプロファイルを求める
    !

    
    !暗黙の型宣言禁止
    implicit none
    
    real(DP), intent(in) :: a_MolFrIni(1:ncmax)
    real(DP), intent(in) :: Humidity
    real(DP), intent(in) :: z_Temp(kmin:kmax)
    real(DP), intent(in) :: z_Press(kmin:kmax)
    real(DP), intent(out):: za_MolFr(kmin:kmax, 1:ncmax)
    
    real(DP)             :: DelMolFr
    integer             :: k, s
    

    !-----------------------------------------------------------
    ! 配列の初期化
    !-----------------------------------------------------------
    do s = 1, ncmax
      za_MolFr(:,s) = a_MolFrIni(s) 
    end do

    !-----------------------------------------------------------
    ! 断熱減率 dT/dz の計算. 
    !-----------------------------------------------------------
    do k = 1, kmax

      za_MolFr(k,:) = za_MolFr(k-1,:)
      
      !------------------------------------------------------------
      !NH4SH 以外の化学種の平衡条件
      !------------------------------------------------------------
      do s = 1, CondNum

        !モル比を求める
        !モル比は前のステップでのモル比を超えることはない
        za_MolFr(k,IdxCG(s)) =                                 &
          & min(                                                &
          &       za_MolFr(k-1,IdxCG(s)),                      &
          &       SvapPress( SpcWetID(IdxCC(s)), z_Temp(k) ) &
          &        * Humidity / z_Press(k)                      &
          &      )
        
      end do

      !------------------------------------------------------------
      !NH4SH の平衡条件
      !------------------------------------------------------------
      if ( IdxNH3 /= 0 ) then 
        
        !モル比の変化. 
        !とりあえず NH4SH に対する飽和比は 1.0 とする(手抜き...).
        DelMolFr =                                            &
          & max (                                             &
          &    DelMolFrNH4SH(                                 &
          &      z_Temp(k), z_Press(k),                       &
          &      za_MolFr(k,IdxNH3), za_MolFr(k,IdxH2S), Humidity   &
          &     ),                                            &
          &    0.0d0                                          &
          &  )
        
        za_MolFr(k,IdxNH3) = za_MolFr(k,IdxNH3) - DelMolFr 
        za_MolFr(k,IdxH2S) = za_MolFr(k,IdxH2S) - DelMolFr
      end if
      
    end do
  end subroutine ECCM_MolFr

!!!------------------------------------------------------------------------------!!!
  subroutine ECCM_DTempDZ( Temp, Press, DelZ, MolFr, DTempDZ )
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(DP), intent(in) :: Temp
    real(DP), intent(in) :: Press
    real(DP), intent(in) :: DelZ
    real(DP), intent(inout) :: MolFr(0:ncmax)    !モル分率
    real(DP), intent(out):: DTempDZ
    real(DP)            :: ReactHeat
    real(DP)            :: Heat(ncmax)
    real(DP)            :: DelMolFr
    real(DP)            :: SatPress
    real(DP)            :: VapPress
    real(DP)            :: Humidity
    real(DP)            :: A, B
    integer             :: s

    !初期化
    DTempDZ      = 0.0d0
    ReactHeat    = 0.0d0
    Heat         = 0.0d0
    DelMolFr     = 0.0d0
    SatPress     = 0.0d0
    VapPress     = 0.0d0

    !------------------------------------------------------------
    !NH4SH 以外の化学種の平衡条件
    !------------------------------------------------------------
    do s = 1, CondNum      
      
      !飽和蒸気圧
      SatPress = SvapPress( SpcWetID(IdxCC(s)), Temp )
      
      !潜熱. 
      Heat(IdxCG(s)) = LatentHeatPerMol( SpcWetID(IdxCC(s)), Temp )
      
      !元々のモル分率を用いて現在の蒸気圧を計算
      VapPress = MolFr(IdxCG(s)) * Press
      
      !飽和蒸気圧から凝結の有無を決める
      if ( VapPress < SatPress ) then         
        !凝結していないので潜熱なし.
        Heat(IdxCG(s)) = 0.0d0          

      else      

        !飽和蒸気圧と圧力から現在のモル比を計算
        MolFr(IdxCG(s)) = max(SatPress / Press, 1.0d-16)

      end if
    end do
    
    !------------------------------------------------------------
    !NH4SH の平衡条件
    !------------------------------------------------------------
    if ( IdxNH3 /= 0 ) then 
      
      Humidity = 1.0d0
      DelMolFr =                                            &
        & max (                                             &
        &    DelMolFrNH4SH(                                 &
        &         Temp, Press, MolFr(IdxNH3), MolFr(IdxH2S),&
        &         Humidity                                  &
        &      ),                                           &
        &    0.0d0                                          &
        &  )
      MolFr(IdxNH3) = MolFr(IdxNH3) - DelMolFr
      MolFr(IdxH2S) = MolFr(IdxH2S) - DelMolFr

      ReactHeat = ReactHeatNH4SHPerMol * DelMolFr
    end if
    
    !------------------------------------------------------------
    !温度勾配を計算
    !------------------------------------------------------------
    !係数. 温度 Temp(i) で評価
    A = dot_product( Heat(1:ncmax), MolFr(1:ncmax)) &
      &  / ( GasRUniv * Temp )
    B = dot_product(( Heat(1:ncmax) ** 2.0d0), MolFr(1:ncmax)) &
      &  / ( CpDryMol * GasRUniv * ( Temp ** 2.0d0 ) )
    
    !断熱温度減率
    DTempDZ = - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B))  &
      &       + ReactHeat / (CpDryMol * DelZ)
    
  end subroutine ECCM_DTempDZ

!!!------------------------------------------------------------------------------!!!
  subroutine ECCM_Stab( xyz_PTemp, xyz_Exner, xyzf_QMix )

    use gridset,only: imin,       &! 配列の X 方向の下限
      &                 imax,       &! 配列の X 方向の上限 
      &                 jmin,       &! 配列の Y 方向の下限
      &                 jmax,       &! 配列の Y 方向の上限 
      &                 kmin,       &! 配列の Z 方向の下限
      &                 kmax,       &! 配列の Z 方向の上限 
      &                 ncmax          !
    use basicset,only:  &
      &                 xyz_ExnerBZ,   &!
      &                 xyz_PTempBZ, &! 
      &                 xyzf_QMixBZ
    use constants,only: MolWtDry,      &!
      &                 CpDry,         &!
      &                 Grav          
    use composition,only:  MolWtWet
    use axesset, only : xyz_avr_xyr
    use xyz_deriv_module,only : xyr_dz_xyz
    
    implicit none

    real(DP), intent(in)  :: xyz_PTemp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyz_Exner(imin:imax,jmin:jmax,kmin:kmax)
    real(DP), intent(in)  :: xyzf_QMix(imin:imax,jmin:jmax,kmin:kmax, ncmax)

    real(DP) :: xyz_Stab(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_StabTemp(imin:imax,jmin:jmax,kmin:kmax)
    real(DP) :: xyz_StabMolWt(imin:imax,jmin:jmax,kmin:kmax)

    real(DP)   :: xyzf_MolFrAll(imin:imax,jmin:jmax,kmin:kmax,ncmax)
    real(DP)   :: xyz_TempAll(imin:imax,jmin:jmax,kmin:kmax)
    real(DP)   :: xyz_MolWtWet(imin:imax,jmin:jmax,kmin:kmax)
    integer    :: i, j, k, s

    xyz_TempAll = (xyz_PTemp + xyz_PTempBZ) * (xyz_Exner + xyz_ExnerBZ)
    do s = 1, ncmax
      xyzf_MolFrAll(:,:,:,s) =                          &
        &   (xyzf_QMix(:,:,:,s) + xyzf_QMixBZ(:,:,:,s)) &
        &   * MolWtDry / MolWtWet(s) 
    end do
    
    do k = kmin, kmax
    do j = jmin, jmax
    do i = imin, imax
      xyz_MolWtWet(i,j,k) = &
        &     dot_product( MolWtWet(1:GasNum), xyzf_MolFrAll(i,j,k,1:GasNum) )
    end do
    end do
    end do
    
    xyz_StabTemp =                                               &
      &         Grav / xyz_TempAll                               &
      &           * (   xyz_avr_xyr( xyr_dz_xyz( xyz_TempAll ) ) &
      &               + Grav / CpDry ) 
    xyz_StabMolWt =                                              &
      &       - Grav * xyz_avr_xyr( xyr_dz_xyz( xyz_MolWtWet ) ) &
      &         / MolWtDry 
    xyz_Stab = xyz_StabTemp + xyz_StabMolWt

    where (xyz_Stab < 1.0d-7) 
      xyz_Stab = 1.0d-7
    end where

  end subroutine ECCM_Stab

end module ECCM