!= Module MoistAdjust
!
! Authors::   SUGIYAMA Ko-ichiro
! Version::   $Id: moistadjust.f90,v 1.14 2006/09/30 05:53:54 odakker Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!湿潤飽和調節法
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

Module MoistAdjust
  !
  !湿潤飽和調節法を行うためのパッケージ型モジュール
  !  * 飽和蒸気圧を使う化学種は, 一元的に実行する
  !  * 化学反応については, それぞれの化学反応毎に実行する. 

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

  use gridset, only:  SpcNum,           &!化学種の数
    &                 DimXMin,          &! x 方向の配列の下限
    &                 DimXMax,          &! x 方向の配列の上限
    &                 DimZMin,          &! z 方向の配列の下限
    &                 DimZMax            ! z 方向の配列の上限
  use basicset, only: CpDry,            &
    &                 MolWtWet,         &
    &                 SpcWetID,         &
    &                 SpcWetSymbol,     &
    &                 xz_ExnerBasicZ,   &
    &                 xz_PotTempBasicZ, &
    &                 xza_MixRtBasicZ,  &  
    &                 PressBasis,       &!温位の基準圧力         [Pa]
    &                 CpDry,            &!乾燥成分の平均定圧比熱 [J/K kg]
    &                 MolWtDry,         &!乾燥成分の平均分子量   [kg/mol]
    &                 GasRDry            !乾燥成分の気体定数     [J/K kg]
  use average,  only: xz_avr_pz, xz_avr_xr, pz_avr_xz, xr_avr_xz
  use ChemData, only: ChemData_OneSpcID
  use ChemCalc, only: xz_SvapPress, xz_LatentHeat, ReactHeatNH4SH
  use MoistFunc,only: xz_DMixRtSatDPotTemp, xz_DelMixRtNH4SH
  use StoreSet, only: StoreCond

  !暗黙の型宣言禁止
  implicit none
  
  !属性の指定
  private

  !関数の属性を public に変更
  public MoistAdjust_Init        !初期化ルーチン
  public MoistAdjustPrm          !変数を出力. 解析ルーチンで利用するため
  public MoistAdjustSvapPress    !飽和蒸気圧を用いた飽和湿潤調節(簡易版)
  public MoistAdjustNH4SH        !化学反応の圧平衡定数を用いた飽和湿潤調節

  !変数定義
  integer      :: LoopNum      = 0
  integer      :: GasNum(10)   = 0
  integer      :: CloudNum(10) = 0
  integer      :: NH3Num   = 0
  integer      :: H2SNum   = 0
  integer      :: NH4SHNum = 0
  
  save LoopNum, CloudNum, GasNum
  save NH3Num, H2SNum, NH4SHNum
  
contains

!!!------------------------------------------------------------------!!!
  subroutine MoistAdjust_Init( )
    !
    !計算に利用する化学種の ID の対を取得
    !  * 気相と凝縮相(雲粒)の ID の対. 
    !  * NH4SH 生成反応を計算するための, NH4SH と NH3, H2S の ID.

    !暗黙の型宣言禁止
    implicit none

    !変数定義
    integer                  :: s
    integer                  :: n1

    !-----------------------------------------------------------
    ! 雲粒と気体の ID の組を作る
    !-----------------------------------------------------------
    !初期化
    LoopNum = 0

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

      ! NH4SH については無視
      if ( trim(SpcWetSymbol(s)) == 'NH4SH-s-Cloud' ) then 
        cycle SelectCloud
      end if

      !'Cloud' という文字列が含まれるものの個数を数える
      n1 = index(SpcWetSymbol(s), '-Cloud' )
      if (n1 /= 0) then
        LoopNum          = LoopNum + 1
        CloudNum(LoopNum)= s
        GasNum(LoopNum)  = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID(SpcWetSymbol(s)(1:n1-3) // '-g'))
      end if
    end do SelectCloud
    
    !-----------------------------------------------------------
    ! 硫化アンモニウム, およびアンモニアと硫化水素の ID を取得
    !-----------------------------------------------------------
    NH3Num   = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH3-g'))
    H2SNum   = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('H2S-g'))
    NH4SHNum = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH4SH-s'))

    !-----------------------------------------------------------
    ! 確認
    !-----------------------------------------------------------
    if ( LoopNum == 0 ) then 
      call MessageNotify( "W", &
        & "MoistAdjust: ", "CloudNum = 0, please comment out of MoistAdjust" )
!      stop
    end if

    call MessageNotify( "M", "MoistAdjust_Init", "LoopNum=%d", i=(/LoopNum/) )
    call MessageNotify( "M", "MoistAdjust_Init", "CloudNum=%d", i=(/CloudNum/) )
    call MessageNotify( "M", "MoistAdjust_Init", "GasNum=%d", i=(/GasNum/) )
    call MessageNotify( "M", "MoistAdjust_Init", "NH3Num=%d", i=(/NH3Num/) )
    call MessageNotify( "M", "MoistAdjust_Init", "H2SNum=%d", i=(/H2SNum/) )
    call MessageNotify( "M", "MoistAdjust_Init", "NH4SHNum=%d", i=(/NH4SHNum/) )

  end subroutine MoistAdjust_Init


  subroutine MoistAdjustPrm( M1, M2, M3, M4, M5, M6 )
    
    implicit none
    integer, intent(out) ::  M1, M2(10), M3(10), M4, M5, M6

    M1 = LoopNum
    M2 = CloudNum
    M3 = GasNum    
    M4 = NH3Num
    M5 = H2SNum
    M6 = NH4SHNum

!    write(*,*) M1, M2, M3, M4, M5, M6 
  end subroutine MoistAdjustPrm


!!!------------------------------------------------------------------!!!
  subroutine MoistAdjustSvapPress(xz_Exner, xz_PotTemp, xza_MixRt)
    !
    ! 飽和蒸気圧を用いた湿潤飽和調整法の実行
    ! この副プログラムでは, 予め決めた回数だけ反復改良を行う. 
    !

    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8),intent(in)   :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !エクスナー関数
    real(8),intent(inout):: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !温位
    real(8),intent(inout):: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                          !混合比
    integer, parameter   :: ItrNum = 4                    !反復改良の回数

    real(8) :: xz_MixRtV_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_MixRtV_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_MixRtC_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_MixRtC_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_PotTemp_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_PotTemp_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_MixRtSat(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_Cond(DimXMin:DimXMax, DimZMin:DimZMax)  
    real(8) :: xz_Evap(DimXMin:DimXMax, DimZMin:DimZMax)  
    real(8) :: xz_Gamma(DimXMin:DimXMax, DimZMin:DimZMax) 
    integer :: i, s                                       ! 添え字  

    !---------------------------------------------------------------------
    ! 初期化
    !---------------------------------------------------------------------
    xz_MixRtV_pre  = 0.0d0
    xz_MixRtV_nxt  = 0.0d0
    xz_MixRtC_pre  = 0.0d0
    xz_MixRtC_nxt  = 0.0d0
    xz_PotTemp_pre = 0.0d0
    xz_PotTemp_nxt = 0.0d0
    
    xz_ExnerAll    = 0.0d0
    xz_TempAll     = 0.0d0
    xz_DelMixRt    = 0.0d0
    xz_MixRtSat    = 0.0d0 
    xz_Gamma       = 0.0d0
    
    !---------------------------------------------------------------------
    ! 湿潤飽和調節法の実行
    !   ループを回すのは, 雲についてだけ.  
    !---------------------------------------------------------------------
    LoopSvapPress: do s = 1, LoopNum
      
      !湿潤飽和法では圧力は変化させない. 
      xz_ExnerAll = xz_Exner + xz_ExnerBasicZ
          
      !今までに得られた擾乱成分の値を暫定量とみなす. 添え字を追加
      xz_MixRtV_pre  = xza_MixRt(:,:,GasNum(s))   + xza_MixRtBasicZ(:,:,GasNum(s))
      xz_MixRtC_pre  = xza_MixRt(:,:,CloudNum(s)) + xza_MixRtBasicZ(:,:,CloudNum(s))
      xz_PotTemp_pre = xz_PotTemp
          
      Adjusting: do i = 1, ItrNum
        !---------------------------------------------------------------
        ! 飽和蒸気圧から飽和混合比を求める
        !---------------------------------------------------------------
        !温度
        xz_TempAll = ( xz_PotTemp_pre + xz_PotTempBasicZ ) * xz_ExnerAll

        !飽和蒸気圧から飽和混合比を計算(基本場からの差). 
        xz_MixRtSat =                                                    &
          & xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll)              &
          &  * MolWtWet(CloudNum(s))                                     &
          &  / (MolWtDry * PressBasis * (xz_ExnerAll ** (CpDry / GasRDry)) )

!        write(*,*) "A", minval( xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll) ), maxval( xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll) )
!        write(*,*) "B", minval( xz_MixRtSat ), maxval( xz_MixRtSat )
        
        !規格化された潜熱
        xz_Gamma = xz_LatentHeat(SpcWetID(CloudNum(s)), xz_TempAll) &
          &        / (xz_ExnerAll * CpDry)
        
!        write(*,*) "C", minval( xz_Gamma ), maxval( xz_Gamma )

        !凝結量を求める. 
        !  凝結が生じる場合には, xz_MixRtV_pre - xz_MixRtSat は必ず正の値となる.
        !  蒸発が生じる場合には, 蒸発量は - MixRtC を超えることはない. 
        xz_DelMixRt = &
          & ( xz_MixRtV_pre - xz_MixRtSat ) &
          &   / (1.0d0 + xz_Gamma * xz_DMixRtSatDPotTemp( &
          &           SpcWetID(CloudNum(s)), MolWtWet(CloudNum(s)), xz_TempAll, xz_ExnerAll & 
          &        ) )

!        Write(*,*) "D1", &
!          & minval( xz_DMixRtSatDPotTemp( &
!          &           SpcWetID(CloudNum(s)), MolWtWet(CloudNum(s)), xz_TempAll, xz_ExnerAll & 
!          &        ) ),      &
!          & maxval( xz_DMixRtSatDPotTemp( &
!          &           SpcWetID(CloudNum(s)), MolWtWet(CloudNum(s)), xz_TempAll, xz_ExnerAll & 
!          &        ) )
!        write(*,*) "D2", minval( xz_MixRtV_pre - xz_MixRtSat ), maxval( xz_MixRtV_pre - xz_MixRtSat )
!        write(*,*) "D3", minval( xz_DelMixRt ), maxval( xz_DelMixRt )

        xz_Cond = max( 0.0d0, min( xz_MixRtV_pre,   xz_DelMixRt ) )
        xz_Evap = max( 0.0d0, min( xz_MixRtC_pre, - xz_DelMixRt ) )

!        write(*,*) "E", minval( xz_Cond ), maxval( xz_Cond )
!        write(*,*) "F", minval( xz_Evap ), maxval( xz_Evap )

        !より真に近い値を計算
        xz_PotTemp_nxt = xz_PotTemp_pre + xz_Gamma * ( xz_Cond - xz_Evap )
        xz_MixRtV_nxt  = xz_MixRtV_pre  - xz_Cond + xz_Evap
        xz_MixRtC_nxt  = xz_MixRtC_pre  + xz_Cond - xz_Evap
        
        !繰り返しのための変数定義
        xz_PotTemp_pre = xz_PotTemp_nxt
        xz_MixRtV_pre  = xz_MixRtV_nxt
        xz_MixRtC_pre  = xz_MixRtC_nxt

      end do Adjusting

      call StoreCond( xz_PotTemp_nxt - xz_PotTemp ) 
      
      xz_PotTemp                 = xz_PotTemp_nxt                 
      xza_MixRt(:,:,GasNum(s))   = xz_MixRtV_nxt - xza_MixRtBasicZ(:,:,GasNum(s))
      xza_MixRt(:,:,CloudNum(s)) = xz_MixRtC_nxt - xza_MixRtBasicZ(:,:,CloudNum(s))
      
    end do LoopSvapPress
    
  end subroutine MoistAdjustSvapPress


!!!--------------------------------------------------------------------------!!!
  subroutine MoistAdjustNH4SH(xz_Exner, xz_PotTemp, xza_MixRt ) 
    !
    ! NH3 + H2S --> NH4SH の生成反応の圧平衡定数 Kp を用いた飽和湿潤調節法
    !
    
    !暗黙の型宣言禁止
    implicit none
    
    !変数定義
    real(8),intent(in)   :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !エクスナー関数
    real(8),intent(inout):: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !温位
    real(8),intent(inout):: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                          !凝縮成分の混合比
    real(8) :: xz_PotTemp_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_PotTemp_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_MixRtNH3_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_MixRtNH3_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_MixRtH2S_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_MixRtH2S_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_MixRtNH4SH_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_MixRtNH4SH_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    
    real(8) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_Gamma(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_Cond(DimXMin:DimXMax, DimZMin:DimZMax)  
    real(8) :: xz_Evap(DimXMin:DimXMax, DimZMin:DimZMax)  
    
    integer            :: i
    integer, parameter :: ItrNum = 2
            
    !---------------------------------------------------------------------
    ! 初期化
    !---------------------------------------------------------------------
    xz_PotTemp_nxt = 0.0d0
    xz_MixRtNH3_nxt   = 0.0d0
    xz_MixRtH2S_nxt   = 0.0d0
    xz_MixRtNH4SH_nxt = 0.0d0
    
    xz_Gamma     = 0.0d0
    xz_DelMixRt  = 0.0d0

    !湿潤飽和法では圧力は変化させない. 
    xz_ExnerAll = xz_Exner + xz_ExnerBasicZ
    xz_PressAll = PressBasis * (xz_ExnerAll ** (CpDry / GasRDry))
        
    !今までに得られた擾乱成分の値を暫定量とみなす. 添え字を追加
    xz_MixRtNH3_pre(:,:)    = xza_MixRt(:,:,NH3Num)   + xza_MixRtBasicZ(:,:,NH3Num)
    xz_MixRtH2S_pre(:,:)    = xza_MixRt(:,:,H2SNum)   + xza_MixRtBasicZ(:,:,H2SNum)
    xz_MixRtNH4SH_pre(:,:)  = xza_MixRt(:,:,NH4SHNum) + xza_MixRtBasicZ(:,:,NH4SHNum)
    xz_PotTemp_pre          = xz_PotTemp
    
    AdjustNH4SH: do i = 1, ItrNum
      !---------------------------------------------------------------
      ! 変数の初期化
      !---------------------------------------------------------------
      !温度
      xz_TempAll = ( xz_PotTemp_pre + xz_PotTempBasicZ ) * xz_ExnerAll
      
      !規格化された反応熱 (NH4SH １kg に対する熱量)
      xz_Gamma = ReactHeatNH4SH / ( xz_ExnerAll * CpDry )

      !NH4SH の生成量
      xz_DelMixRt =                                                           &
        &    xz_DelMixRtNH4SH(                                                &
        &          xz_TempAll, xz_PressAll, xz_MixRtNH3_pre, xz_MixRtH2S_pre, &
        &          MolWtWet(NH3Num), MolWtWet(H2SNum)                         &
        &      )

      xz_Cond = max( 0.0d0, xz_DelMixRt )
      xz_Evap = max( 0.0d0, min( - xz_DelMixRt, xz_MixRtNH4SH_pre ) )

      !---------------------------------------------------------------
      ! より真に近い値を求める飽和蒸気圧から飽和混合比を求める
      !---------------------------------------------------------------
      ! NH4SH の混合比を修正
      xz_MixRtNH4SH_nxt  = xz_MixRtNH4SH_pre + xz_Cond - xz_Evap
      
      ! DelPress を元に, NH3 と H2S の混合比を修正
      xz_MixRtNH3_nxt = xz_MixRtNH3_pre - ( xz_Cond - xz_Evap )      &
        &                 * MolWtWet(NH3Num) / MolWtWet(NH4SHNum)
      xz_MixRtH2S_nxt = xz_MixRtH2S_pre - ( xz_Cond - xz_Evap )      &
        &                 * MolWtWet(H2SNum) / MolWtWet(NH4SHNum)
          
      !温位を修正
      xz_PotTemp_nxt = xz_PotTemp_pre + xz_Gamma * ( xz_Cond - xz_Evap )
      
      !ループを回すための変数変化
      xz_PotTemp_pre     = xz_PotTemp_nxt
      xz_MixRtNH3_pre   = xz_MixRtNH3_nxt 
      xz_MixRtH2S_pre   = xz_MixRtH2S_nxt 
      xz_MixRtNH4SH_pre = xz_MixRtNH4SH_nxt 

    end do AdjustNH4SH

    call StoreCond( xz_PotTemp_nxt - xz_PotTemp ) 

    xz_PotTemp              = xz_PotTemp_nxt                 
    xza_MixRt(:,:,NH3Num)   = xz_MixRtNH3_nxt   - xza_MixRtBasicZ(:,:,NH3Num)
    xza_MixRt(:,:,H2SNum)   = xz_MixRtH2S_nxt   - xza_MixRtBasicZ(:,:,H2SNum)
    xza_MixRt(:,:,NH4SHNum) = xz_MixRtNH4SH_nxt - xza_MixRtBasicZ(:,:,NH4SHNum)
    
  end subroutine MoistAdjustNH4SH

  
end Module MoistAdjust
