!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Module MoistAdjust
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!湿潤飽和調節法
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

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

  !モジュール読み込み
  use gridset, only:  SpcNum,           &!化学種の数
    &                 DimXMin,          &! x 方向の配列の下限
    &                 DimXMax,          &! x 方向の配列の上限
    &                 DimZMin,          &! z 方向の配列の下限
    &                 DimZMax,          &! z 方向の配列の上限
    &                 RegXMin,          &! x 方向の物理領域の下限
    &                 RegXMax,          &! x 方向の物理領域の上限
    &                 RegZMin,          &! z 方向の物理領域の下限
    &                 RegZMax            ! z 方向の物理領域の上限
  use basicset, only: CpDry,            &
    &                 MolWtWet,         &
    &                 SpcWetID,         &
    &                 SpcWetSymbol,     &
    &                 xz_ExnerBasicZ,   &
    &                 xz_PotTempBasicZ, &
    &                 xz_MixRtBasicZ  
  use average,  only: xz_avr_pz, xz_avr_xr, pz_avr_xz, xr_avr_xz
  use ChemData, only: ChemData_OneSpcID
  use ChemCalc, only: SvapPress, LatentHeatPerMass,  &
    &                 xz_SvapPress, xz_LatentHeatPerMass, &
    &                 ReactHeatNH4SHPerMass
  use MoistFunc,only: Vap2MixRt, DMixRtSatDPotTemp, MixRt2Vap, &
    &                 xz_Vap2MixRt, xz_DMixRtSatDPotTemp, xz_MixRt2Vap, &
    &                 Exner2Press, DelMixRtNH4SH, &
    &                 xz_Exner2Press, xz_DelMixRtNH4SH

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

  !関数の属性を public に変更
  public MoistAdjust_Init        !初期化ルーチン
  public MoistAdjustSvapPress    !飽和蒸気圧を用いた飽和湿潤調節
  public MoistAdjustSvapPress2   !飽和蒸気圧を用いた飽和湿潤調節(簡易版)
  public MoistAdjustNH4SH        !化学反応の圧平衡定数を用いた飽和湿潤調節
  public MoistAdjustNH4SH2       !化学反応の圧平衡定数を用いた飽和湿潤調節

  !変数定義
  real(8)      :: MinError = 1.0d-12
  integer      :: LoopNum      = 0
  integer      :: GasNum(10)   = 0
  integer      :: CloudNum(10) = 0
  integer      :: NH3Num   = 0
  integer      :: H2SNum   = 0
  integer      :: NH4SHNum = 0
  
  save MinError
  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 
      write(*,*) "MoistAdjust: CloudNum = 0, please comment out of MoistAdjust"
!      stop
    end if

    write(*,*) "MoistAdjust_Init, LoopNum:  ", LoopNum
    write(*,*) "MoistAdjust_Init, CloudNum: ", CloudNum
    write(*,*) "MoistAdjust_Init, GasNum:   ", GasNum    
    write(*,*) "MoistAdjust_Init, NH3Num:   ", NH3Num
    write(*,*) "MoistAdjust_Init, H2SNum:   ", H2SNum
    write(*,*) "MoistAdjust_Init, NH4SHNum: ", NH4SHNum

  end subroutine MoistAdjust_Init


!!!------------------------------------------------------------------!!!
  subroutine MoistAdjustSvapPress(xz_Exner, xz_PotTemp, xz_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):: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                          !混合比
    real(8) :: MixRtV_pre
    real(8) :: MixRtV_nxt
    real(8) :: MixRtC_pre
    real(8) :: MixRtC_nxt
    real(8) :: PotTemp_pre
    real(8) :: PotTemp_nxt    
    real(8) :: ExnerAll
    real(8) :: TempAll
    real(8) :: DelMixRt
    real(8) :: MixRtSat
    real(8) :: Gamma

    integer :: i, k, s, num          ! 添え字  
    
    !---------------------------------------------------------------------
    ! 初期化
    !---------------------------------------------------------------------
    MixRtV_pre  = 0.0d0
    MixRtV_nxt  = 0.0d0
    MixRtC_pre  = 0.0d0
    MixRtC_nxt  = 0.0d0
    PotTemp_pre = 0.0d0
    PotTemp_nxt = 0.0d0
    
    ExnerAll    = 0.0d0
    TempAll     = 0.0d0
    DelMixRt    = 0.0d0
    MixRtSat    = 0.0d0 
    Gamma       = 0.0d0
    
    
    num   = 0
    
    !ループを回すのは, 雲についてだけ.  
    LoopSvapPress: do s = 1, LoopNum
!      write(*,*) 'pre: ', minval(xz_MixRt(:,:,GasNum(s))), maxval(xz_MixRt(:,:,GasNum(s)))
!      write(*,*) 'pre: ', minval(xz_MixRt(:,:,CloudNum(s))), maxval(xz_MixRt(:,:,CloudNum(s)))
            
      do k = RegZMin, RegZMax
        do i = RegXMin, RegXMax
          
          !湿潤飽和法では圧力は変化させない. 
          ExnerAll = xz_Exner(i,k) + xz_ExnerBasicZ(i,k)
          
          !今までに得られた擾乱成分の値を暫定量とみなす. 添え字を追加
          MixRtV_pre  = xz_MixRt(i,k,GasNum(s))   + xz_MixRtBasicZ(i,k,GasNum(s))
          MixRtC_pre  = xz_MixRt(i,k,CloudNum(s)) + xz_MixRtBasicZ(i,k,CloudNum(s))
          PotTemp_pre = xz_PotTemp(i,k)
          
          !ループ回数の初期化
          num = 0
          
          Adjusting: do 
            !---------------------------------------------------------------
            ! 飽和蒸気圧から飽和混合比を求める
            !---------------------------------------------------------------
            !温度
            TempAll = (PotTemp_pre + xz_PotTempBasicZ(i,k)) * ExnerAll
            
            !飽和蒸気圧から飽和混合比を計算(基本場からの差). 
            MixRtSat = Vap2MixRt(MolWtWet(CloudNum(s)), SvapPress(SpcWetID(CloudNum(s)), TempAll), ExnerAll)
            
            !規格化された潜熱
            Gamma = LatentHeatPerMass(SpcWetID(CloudNum(s)), TempAll) / (ExnerAll * CpDry)
            
            !---------------------------------------------------------------
            ! 雲粒への変換を行う
            ! 水蒸気分圧が飽和蒸気圧より多い or 雲粒混合比が正
            !---------------------------------------------------------------
            if (MixRtV_pre - MixRtSat > MinError .OR. MixRtC_pre > MinError) then
              
              !凝結量を求める
              DelMixRt = max(- MixRtC_pre, ( MixRtV_pre - MixRtSat ) )                          &
                & / (1.0d0 + Gamma * DMixRtSatDPotTemp(SpcWetID(CloudNum(s)), MolWtWet(CloudNum(s)), TempAll, ExnerAll))
              
              PotTemp_nxt = PotTemp_pre + Gamma * DelMixRt
              MixRtV_nxt  = MixRtV_pre - DelMixRt
              MixRtC_nxt  = MixRtC_pre + DelMixRt
              
              !温位が収束していたら繰り返しを中止する
              if (abs((PotTemp_nxt - PotTemp_pre) / PotTemp_nxt) <= MinError ) then 
                xz_PotTemp(i,k)          = PotTemp_nxt                 
                xz_MixRt(i,k,GasNum(s))   = MixRtV_nxt - xz_MixRtBasicZ(i,k,GasNum(s))
                xz_MixRt(i,k,CloudNum(s)) = MixRtC_nxt - xz_MixRtBasicZ(i,k,CloudNum(s))
                exit Adjusting
              end if
            
              !ループ回数が回りすぎたらエラーを吐いて続ける
              if ( num > 20 ) then 
                write(*,*) "Looping is too much, num > 20"
                xz_PotTemp(i,k)          = PotTemp_nxt     
                xz_MixRt(i,k,GasNum(s))   = MixRtV_nxt - xz_MixRtBasicZ(i,k,GasNum(s))
                xz_MixRt(i,k,CloudNum(s)) = MixRtC_nxt - xz_MixRtBasicZ(i,k,CloudNum(s))
                exit Adjusting
              end if
              
            else 
              
              !収束したとして終了
              !雲粒混合比が負になったときは, 雲粒混合比を正にする. 
              xz_PotTemp(i,k)          = PotTemp_pre - Gamma * MixRtC_pre
              xz_MixRt(i,k,GasNum(s))   = MixRtV_pre + MixRtC_pre - xz_MixRtBasicZ(i,k,GasNum(s))
              xz_MixRt(i,k,CloudNum(s)) = 0.0d0          
              exit Adjusting
            end if
            
            !---------------------------------------------------------------
            ! 繰り返しのための変数定義
            !---------------------------------------------------------------
            PotTemp_pre = PotTemp_nxt
            MixRtV_pre  = MixRtV_nxt
            MixRtC_pre  = MixRtC_nxt
            num         = num + 1
            
          end do Adjusting
          
        end do
      end do
 
!      write(*,*) 'nxt: ', minval(xz_MixRt(:,:,GasNum(s))), maxval(xz_MixRt(:,:,GasNum(s)))
!      write(*,*) 'nxt: ', minval(xz_MixRt(:,:,CloudNum(s))), maxval(xz_MixRt(:,:,CloudNum(s)))
    end do LoopSvapPress

    
  end subroutine MoistAdjustSvapPress

  
!!!------------------------------------------------------------------!!!
  subroutine MoistAdjustSvapPress2(xz_Exner, xz_PotTemp, xz_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):: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                          !混合比
    integer, parameter   :: ItrNum = 3                    !反復改良の回数
    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_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  = xz_MixRt(:,:,GasNum(s))   + xz_MixRtBasicZ(:,:,GasNum(s))
      xz_MixRtC_pre  = xz_MixRt(:,:,CloudNum(s)) + xz_MixRtBasicZ(:,:,CloudNum(s))
      xz_PotTemp_pre = xz_PotTemp
          
      Adjusting: do i = 0, ItrNum
        !---------------------------------------------------------------
        ! 飽和蒸気圧から飽和混合比を求める
        !---------------------------------------------------------------
        !温度
        xz_TempAll = (xz_PotTemp_pre + xz_PotTempBasicZ) * xz_ExnerAll

        !飽和蒸気圧から飽和混合比を計算(基本場からの差). 
        xz_MixRtSat = xz_Vap2MixRt( MolWtWet(CloudNum(s)), xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll), xz_ExnerAll )
        
        !規格化された潜熱
        xz_Gamma = xz_LatentHeatPerMass(SpcWetID(CloudNum(s)), xz_TempAll) / (xz_ExnerAll * CpDry)
        
        !凝結量を求める. 
        !  凝結が生じる場合には, xz_MixRtV_pre - xz_MixRtSat は必ず正の値となる.
        !  蒸発が生じる場合には, 蒸発量は - MixRtC を超えることはない. 
        !  merge 関数を用いて蒸発量の上限を定めている. 
        xz_DelMixRt = &
          & merge( &
          &    (xz_MixRtV_pre - xz_MixRtSat), &
          &    - xz_MixRtC_pre, &
          &    (xz_MixRtV_pre - xz_MixRtSat) > - xz_MixRtC_pre &
          &   ) &
          &   / (1.0d0 + xz_Gamma * xz_DMixRtSatDPotTemp( &
          &           SpcWetID(CloudNum(s)), MolWtWet(CloudNum(s)), xz_TempAll, xz_ExnerAll & 
          &        ) )
        
        !より真に近い値を計算
        xz_PotTemp_nxt = xz_PotTemp_pre + xz_Gamma * xz_DelMixRt
        xz_MixRtV_nxt  = xz_MixRtV_pre  - xz_DelMixRt
        xz_MixRtC_nxt  = xz_MixRtC_pre  + xz_DelMixRt
        
        !繰り返しのための変数定義
        xz_PotTemp_pre = xz_PotTemp_nxt
        xz_MixRtV_pre  = xz_MixRtV_nxt
        xz_MixRtC_pre  = xz_MixRtC_nxt
        
      end do Adjusting
      
      !反復改良後の値を真の値として保管
      xz_PotTemp                = xz_PotTemp_nxt                 
      xz_MixRt(:,:,GasNum(s))   = xz_MixRtV_nxt - xz_MixRtBasicZ(:,:,GasNum(s))
      xz_MixRt(:,:,CloudNum(s)) = xz_MixRtC_nxt - xz_MixRtBasicZ(:,:,CloudNum(s))
      
    end do LoopSvapPress
    
  end subroutine MoistAdjustSvapPress2

  
!!!--------------------------------------------------------------------------!!!
  subroutine MoistAdjustNH4SH(xz_Exner, xz_PotTemp, xz_MixRt ) 
    !
    ! 化学反応の圧平衡定数 Kp を用いた飽和湿潤調節法
    ! ただし, 現在は NH3 + H2S --> NH4SH の生成反応にのみ対応
    !
    
    !暗黙の型宣言禁止
    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):: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                          !凝縮成分の混合比
    real(8) :: PotTemp_pre
    real(8) :: PotTemp_nxt
    real(8) :: MixRtNH3_pre
    real(8) :: MixRtNH3_nxt
    real(8) :: MixRtH2S_pre
    real(8) :: MixRtH2S_nxt
    real(8) :: MixRtNH4SH_pre
    real(8) :: MixRtNH4SH_nxt
    
    real(8) :: ExnerAll
    real(8) :: TempAll
    real(8) :: Gamma
    real(8) :: DelMixRt
    
    integer :: i, k, num          ! 添え字  
            
    !---------------------------------------------------------------------
    ! 初期化
    !---------------------------------------------------------------------
    PotTemp_pre = 0.0d0
    PotTemp_nxt = 0.0d0
    MixRtNH3_pre   = 0.0d0
    MixRtNH3_nxt   = 0.0d0
    MixRtH2S_pre   = 0.0d0
    MixRtH2S_nxt   = 0.0d0
    MixRtNH4SH_pre = 0.0d0
    MixRtNH4SH_nxt = 0.0d0
    
    ExnerAll  = 0.0d0
    TempAll   = 0.0d0
    Gamma     = 0.0d0
    DelMixRt  = 0.0d0
    num       = 0
  
    do k = RegZMin, RegZMax
      do i = RegXMin, RegXMax

        !湿潤飽和法では圧力は変化させない. 
        ExnerAll = xz_Exner(i,k) + xz_ExnerBasicZ(i,k)
        
        !今までに得られた擾乱成分の値を暫定量とみなす. 添え字を追加
        MixRtNH3_pre    = xz_MixRt(i,k,NH3Num)   + xz_MixRtBasicZ(i,k,NH3Num)
        MixRtH2S_pre    = xz_MixRt(i,k,H2SNum)   + xz_MixRtBasicZ(i,k,H2SNum)
        MixRtNH4SH_pre  = xz_MixRt(i,k,NH4SHNum) + xz_MixRtBasicZ(i,k,NH4SHNum)
        PotTemp_pre     = xz_PotTemp(i,k)
        
        !ループ回数の初期化
        num = 0
        
        AdjustNH4SH: do
          !---------------------------------------------------------------
          ! 変数の初期化
          !---------------------------------------------------------------
          !温度
          TempAll = (PotTemp_pre + xz_PotTempBasicZ(i,k)) * ExnerAll
          
          !規格化された反応熱 (NH4SH １kg に対する熱量)
          Gamma = ReactHeatNH4SHPerMass / (ExnerAll * CpDry)
          
          DelMixRt = DelMixRtNH4SH(TempAll, Exner2Press(ExnerAll),                         &
            &                      MixRtNH3_pre, MixRtH2S_pre, MixRtNH4SH_pre,             &
            &                      MolWtWet(NH3Num), MolWtWet(H2SNum), MolWtWet(NH4SHNum) )  

          !---------------------------------------------------------------
          ! より真に近い値を求める飽和蒸気圧から飽和混合比を求める
          !---------------------------------------------------------------
          ! NH4SH の混合比を修正
          MixRtNH4SH_nxt  = MixRtNH4SH_pre + DelMixRt
          
          ! DelPress を元に, NH3 と H2S の混合比を修正
          MixRtNH3_nxt = MixRtNH3_pre - DelMixRt * MolWtWet(NH3Num) / MolWtWet(NH4SHNum)
          MixRtH2S_nxt = MixRtH2S_pre - DelMixRt * MolWtWet(H2SNum) / MolWtWet(NH4SHNum)
          
          !温位を修正
          PotTemp_nxt = PotTemp_pre + Gamma * DelMixRt
          
!          if (k > 22 .AND. k < 24 .AND. i == 10) then 
!            write(*,*) 'CHECK'
!            write(*,*) num, i, k, DelMixRt
!            write(*,*) num, i, k, Gamma, Gamma * DelMixRt
!            
!            write(*,*) num, i, k, PotTemp_pre,    " --> ", PotTemp_nxt
!            write(*,*) num, i, k, MixRtNH3_pre,   " --> ", MixRtNH3_nxt
!            write(*,*) num, i, k, MixRtH2S_pre,   " --> ", MixRtH2S_nxt
!            write(*,*) num, i, k, MixRtNH4SH_pre, " --> ", MixRtNH4SH_nxt
!
!            write(*,*) 'pre', MixRtNH3_pre + MixRtH2S_pre + MixRtNH4SH_pre 
!            write(*,*) 'nxt', MixRtNH3_nxt + MixRtH2S_nxt + MixRtNH4SH_nxt 
!          end if
          
          !---------------------------------------------------------------
          ! 収束判定
          !---------------------------------------------------------------
          !温位が収束していたら繰り返しを中止する
          if (abs((PotTemp_nxt - PotTemp_pre) / PotTemp_nxt) <= MinERROR   &
            & .OR. abs(PotTemp_nxt - PotTemp_pre) <= MinERROR             ) then 
            xz_PotTemp(i,k)         = PotTemp_nxt                 
            xz_MixRt(i,k,NH3Num)   = MixRtNH3_nxt - xz_MixRtBasicZ(i,k,NH3Num)
            xz_MixRt(i,k,H2SNum)   = MixRtH2S_nxt - xz_MixRtBasicZ(i,k,H2SNum)
            xz_MixRt(i,k,NH4SHNum) = MixRtNH4SH_nxt - xz_MixRtBasicZ(i,k,NH4SHNum)
            exit AdjustNH4SH
          end if
          
          !ループ回数が回りすぎたらエラーを吐いて続ける
          if ( num > 30 ) then 
            write(*,*) "Looping is too much, num > 30"
            write(*,*) "xz_PotTemp(i,k) ", i, k, PotTemp_nxt, PotTemp_pre
            xz_PotTemp(i,k)         = PotTemp_nxt                 
            xz_MixRt(i,k,NH3Num)   = MixRtNH3_nxt - xz_MixRtBasicZ(i,k,NH3Num)
            xz_MixRt(i,k,H2SNum)   = MixRtH2S_nxt - xz_MixRtBasicZ(i,k,H2SNum)
            xz_MixRt(i,k,NH4SHNum) = MixRtNH4SH_nxt - xz_MixRtBasicZ(i,k,NH4SHNum)
            exit AdjustNH4SH
          end if
          
          !繰り返しのための変数定義
          PotTemp_pre    = PotTemp_nxt
          MixRtNH3_pre   = MixRtNH3_nxt 
          MixRtH2S_pre   = MixRtH2S_nxt 
          MixRtNH4SH_pre = MixRtNH4SH_nxt 
          num            = num + 1
          
        end do AdjustNH4SH
      end do 
    end do 
    
  end subroutine MoistAdjustNH4SH


!!!--------------------------------------------------------------------------!!!
  subroutine MoistAdjustNH4SH2(xz_Exner, xz_PotTemp, xz_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):: xz_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_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_Gamma(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax) 
    
    integer            :: i
    integer, parameter :: ItrNum = 5
            
    !---------------------------------------------------------------------
    ! 初期化
    !---------------------------------------------------------------------
    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_MixRtNH3_pre(:,:)    = xz_MixRt(:,:,NH3Num)   + xz_MixRtBasicZ(:,:,NH3Num)
    xz_MixRtH2S_pre(:,:)    = xz_MixRt(:,:,H2SNum)   + xz_MixRtBasicZ(:,:,H2SNum)
    xz_MixRtNH4SH_pre(:,:)  = xz_MixRt(:,:,NH4SHNum) + xz_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 = ReactHeatNH4SHPerMass / ( xz_ExnerAll * CpDry )

      !NH4SH の生成量
      xz_DelMixRt = xz_DelMixRtNH4SH(xz_TempAll, xz_Exner2Press(xz_ExnerAll),           &
        &                      xz_MixRtNH3_pre, xz_MixRtH2S_pre, xz_MixRtNH4SH_pre,     &
        &                      MolWtWet(NH3Num), MolWtWet(H2SNum), MolWtWet(NH4SHNum) )  

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

    xz_PotTemp             = xz_PotTemp_nxt                 
    xz_MixRt(:,:,NH3Num)   = xz_MixRtNH3_nxt - xz_MixRtBasicZ(:,:,NH3Num)
    xz_MixRt(:,:,H2SNum)   = xz_MixRtH2S_nxt - xz_MixRtBasicZ(:,:,H2SNum)
    xz_MixRt(:,:,NH4SHNum) = xz_MixRtNH4SH_nxt - xz_MixRtBasicZ(:,:,NH4SHNum)
    
  end subroutine MoistAdjustNH4SH2
  
end Module MoistAdjust
