!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Module MoistAdjust
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: moistadjust.f90,v 1.1.1.1 2006/04/25 03:43:58 deepconv Exp $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!$B<>=aK0OBD4@aK!(B
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!== Future Plans
!
!

Module MoistAdjust
  !
  !$B<>=aK0OBD4@aK!$r9T$&$?$a$N%Q%C%1!<%87?%b%8%e!<%k(B
  !  * $BK0OB>x5$05$r;H$&2=3X<o$O(B, $B0l85E*$K<B9T$9$k(B
  !  * $B2=3XH?1~$K$D$$$F$O(B, $B$=$l$>$l$N2=3XH?1~Kh$K<B9T$9$k(B. 

  !$B%b%8%e!<%kFI$_9~$_(B
  use gridset, only:  SpcNum,           &!$B2=3X<o$N?t(B
    &                 DimXMin,          &! x $BJ}8~$NG[Ns$N2<8B(B
    &                 DimXMax,          &! x $BJ}8~$NG[Ns$N>e8B(B
    &                 DimZMin,          &! z $BJ}8~$NG[Ns$N2<8B(B
    &                 DimZMax,          &! z $BJ}8~$NG[Ns$N>e8B(B
    &                 RegXMin,          &! x $BJ}8~$NJ*M}NN0h$N2<8B(B
    &                 RegXMax,          &! x $BJ}8~$NJ*M}NN0h$N>e8B(B
    &                 RegZMin,          &! z $BJ}8~$NJ*M}NN0h$N2<8B(B
    &                 RegZMax            ! z $BJ}8~$NJ*M}NN0h$N>e8B(B
  use basicset, only: CpDry,            &
    &                 MolWtWet,         &
    &                 SpcWetID,         &
    &                 SpcWetSymbol,     &
    &                 xz_ExnerBasicZ,   &
    &                 xz_PotTempBasicZ, &
    &                 xza_MixRtBasicZ,  &  
    &                 PressSfc,         &!$BI8=`05NO(B               [Pa]
    &                 CpDry,            &!$B4%Ag@.J,$NJ?6QDj05HfG.(B [J/K kg]
    &                 MolWtDry,         &!$B4%Ag@.J,$NJ?6QJ,;RNL(B   [kg/mol]
    &                 GasRDry            !$B4%Ag@.J,$N5$BNDj?t(B     [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: SvapPress, LatentHeatPerMass,  &
    &                 xz_SvapPress, xz_LatentHeatPerMass, &
    &                 ReactHeatNH4SHPerMass
  use MoistFunc,only: DMixRtSatDPotTemp, xz_DMixRtSatDPotTemp, &
    &                 DelMixRtNH4SH, xz_DelMixRtNH4SH

  !$B0EL[$N7?@k8@6X;_(B
  implicit none
  
  !$BB0@-$N;XDj(B
  private

  !$B4X?t$NB0@-$r(B public $B$KJQ99(B
  public MoistAdjust_Init        !$B=i4|2=%k!<%A%s(B
  public MoistAdjustSvapPress    !$BK0OB>x5$05$rMQ$$$?K0OB<>=aD4@a(B
  public MoistAdjustSvapPress2   !$BK0OB>x5$05$rMQ$$$?K0OB<>=aD4@a(B($B4J0WHG(B)
  public MoistAdjustNH4SH        !$B2=3XH?1~$N05J?9UDj?t$rMQ$$$?K0OB<>=aD4@a(B
  public MoistAdjustNH4SH2       !$B2=3XH?1~$N05J?9UDj?t$rMQ$$$?K0OB<>=aD4@a(B

  !$BJQ?tDj5A(B
  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( )
    !
    !$B7W;;$KMxMQ$9$k2=3X<o$N(B ID $B$NBP$r<hF@(B
    !  * $B5$Aj$H6E=LAj(B($B1@N3(B)$B$N(B ID $B$NBP(B. 
    !  * NH4SH $B@8@.H?1~$r7W;;$9$k$?$a$N(B, NH4SH $B$H(B NH3, H2S $B$N(B ID.

    !$B0EL[$N7?@k8@6X;_(B
    implicit none

    !$BJQ?tDj5A(B
    integer                  :: s
    integer                  :: n1

    !-----------------------------------------------------------
    ! $B1@N3$H5$BN$N(B ID $B$NAH$r:n$k(B
    !-----------------------------------------------------------
    !$B=i4|2=(B
    LoopNum = 0

    !$B2=3X<o$NCf$+$i1@N3$r:n$k$b$N$rA*$S(B, $B$=$NG[NsE:$(;z$HJ,;RNL$rJ]4I(B.
    SelectCloud: do s = 1, SpcNum

      ! NH4SH $B$K$D$$$F$OL5;k(B
      if ( trim(SpcWetSymbol(s)) == 'NH4SH-s-Cloud' ) then 
        cycle SelectCloud
      end if

      !'Cloud' $B$H$$$&J8;zNs$,4^$^$l$k$b$N$N8D?t$r?t$($k(B
      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
    
    !-----------------------------------------------------------
    ! $BN22=%"%s%b%K%&%`(B, $B$*$h$S%"%s%b%K%"$HN22=?eAG$N(B ID $B$r<hF@(B
    !-----------------------------------------------------------
    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'))

    !-----------------------------------------------------------
    ! $B3NG'(B
    !-----------------------------------------------------------
    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, xza_MixRt ) 
    !
    ! $BK0OB>x5$05$rMQ$$$?<>=aK0OBD4@0K!$N<B9T(B
    !

    !$B0EL[$N7?@k8@6X;_(B
    implicit none
    
    !$BJQ?tDj5A(B
    real(8),intent(in)   :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B%(%/%9%J!<4X?t(B
    real(8),intent(inout):: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B290L(B
    real(8),intent(inout):: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                          !$B:.9gHf(B
    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          ! $BE:$(;z(B  
    
    !---------------------------------------------------------------------
    ! $B=i4|2=(B
    !---------------------------------------------------------------------
    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
    
    !$B%k!<%W$r2s$9$N$O(B, $B1@$K$D$$$F$@$1(B.  
    LoopSvapPress: do s = 1, LoopNum
!      write(*,*) 'pre: ', minval(xza_MixRt(:,:,GasNum(s))), maxval(xza_MixRt(:,:,GasNum(s)))
!      write(*,*) 'pre: ', minval(xza_MixRt(:,:,CloudNum(s))), maxval(xza_MixRt(:,:,CloudNum(s)))
            
      do k = RegZMin, RegZMax
        do i = RegXMin, RegXMax
          
          !$B<>=aK0OBK!$G$O05NO$OJQ2=$5$;$J$$(B. 
          ExnerAll = xz_Exner(i,k) + xz_ExnerBasicZ(i,k)
          
          !$B:#$^$G$KF@$i$l$?>qMp@.J,$NCM$r;CDjNL$H$_$J$9(B. $BE:$(;z$rDI2C(B
          MixRtV_pre  = xza_MixRt(i,k,GasNum(s))   + xza_MixRtBasicZ(i,k,GasNum(s))
          MixRtC_pre  = xza_MixRt(i,k,CloudNum(s)) + xza_MixRtBasicZ(i,k,CloudNum(s))
          PotTemp_pre = xz_PotTemp(i,k)
          
          !$B%k!<%W2s?t$N=i4|2=(B
          num = 0
          
          Adjusting: do 
            !---------------------------------------------------------------
            ! $BK0OB>x5$05$+$iK0OB:.9gHf$r5a$a$k(B
            !---------------------------------------------------------------
            !$B29EY(B
            TempAll = (PotTemp_pre + xz_PotTempBasicZ(i,k)) * ExnerAll
            
            !$BK0OB>x5$05$+$iK0OB:.9gHf$r7W;;(B($B4pK\>l$+$i$N:9(B). 
            MixRtSat =                                                    &
              & SvapPress(SpcWetID(CloudNum(s)), TempAll)                 &
              &  * MolWtWet(CloudNum(s))                                  &
              &  / (MolWtDry * PressSfc * (ExnerAll ** (CpDry / GasRDry)) )
            
            !$B5,3J2=$5$l$?@xG.(B
            Gamma = LatentHeatPerMass(SpcWetID(CloudNum(s)), TempAll) / (ExnerAll * CpDry)
            
            !---------------------------------------------------------------
            ! $B1@N3$X$NJQ49$r9T$&(B
            ! $B?e>x5$J,05$,K0OB>x5$05$h$jB?$$(B or $B1@N3:.9gHf$,@5(B
            !---------------------------------------------------------------
            if (MixRtV_pre - MixRtSat > MinError .OR. MixRtC_pre > MinError) then
              
              !$B6E7kNL$r5a$a$k(B
              DelMixRt = ( 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( DelMixRt < 0.0d0 .AND. i == 50) write(*,*) MixRtV_nxt, MixRtSat, MixRtV_nxt - MixRtSat

              !$B290L$,<}B+$7$F$$$?$i7+$jJV$7$rCf;_$9$k(B
              if (abs((PotTemp_nxt - PotTemp_pre) / PotTemp_nxt) <= MinError ) then 
                xz_PotTemp(i,k)          = PotTemp_nxt                 
                xza_MixRt(i,k,GasNum(s))   = MixRtV_nxt - xza_MixRtBasicZ(i,k,GasNum(s))
                xza_MixRt(i,k,CloudNum(s)) = MixRtC_nxt - xza_MixRtBasicZ(i,k,CloudNum(s))
                exit Adjusting
              end if
            
              !$B%k!<%W2s?t$,2s$j$9$.$?$i%(%i!<$rEG$$$FB3$1$k(B
              if ( num > 20 ) then 
                write(*,*) "Looping is too much, num > 20"
                xz_PotTemp(i,k)          = PotTemp_nxt     
                xza_MixRt(i,k,GasNum(s))   = MixRtV_nxt - xza_MixRtBasicZ(i,k,GasNum(s))
                xza_MixRt(i,k,CloudNum(s)) = MixRtC_nxt - xza_MixRtBasicZ(i,k,CloudNum(s))
                exit Adjusting
              end if
              
            else 
              
              !$B<}B+$7$?$H$7$F=*N;(B
              !$B1@N3:.9gHf$,Ii$K$J$C$?$H$-$O(B, $B1@N3:.9gHf$r@5$K$9$k(B. 
              xz_PotTemp(i,k)          = PotTemp_pre - Gamma * MixRtC_pre
              xza_MixRt(i,k,GasNum(s))   = MixRtV_pre + MixRtC_pre - xza_MixRtBasicZ(i,k,GasNum(s))
              xza_MixRt(i,k,CloudNum(s)) = 0.0d0          
              exit Adjusting
            end if
            
            !---------------------------------------------------------------
            ! $B7+$jJV$7$N$?$a$NJQ?tDj5A(B
            !---------------------------------------------------------------
            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(xza_MixRt(:,:,GasNum(s))), maxval(xza_MixRt(:,:,GasNum(s)))
!      write(*,*) 'nxt: ', minval(xza_MixRt(:,:,CloudNum(s))), maxval(xza_MixRt(:,:,CloudNum(s)))
    end do LoopSvapPress

    
  end subroutine MoistAdjustSvapPress

  
!!!------------------------------------------------------------------!!!
  subroutine MoistAdjustSvapPress2(xz_Exner, xz_PotTemp, xza_MixRt ) 
    !
    ! $BK0OB>x5$05$rMQ$$$?<>=aK0OBD4@0K!$N<B9T(B
    ! $B$3$N%k!<%A%s$G$O(B, $BM=$a7h$a$?2s?t$@$1H?I|2~NI$r9T$&(B. 
    !

    !$B0EL[$N7?@k8@6X;_(B
    implicit none
    
    !$BJQ?tDj5A(B
    real(8),intent(in)   :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B%(%/%9%J!<4X?t(B
    real(8),intent(inout):: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B290L(B
    real(8),intent(inout):: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                          !$B:.9gHf(B
    integer, parameter   :: ItrNum = 5                    !$BH?I|2~NI$N2s?t(B
    real(8) :: xza_MixRtV_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xza_MixRtV_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xza_MixRtC_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xza_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                                       ! $BE:$(;z(B  
    
    !---------------------------------------------------------------------
    ! $B=i4|2=(B
    !---------------------------------------------------------------------
    xza_MixRtV_pre  = 0.0d0
    xza_MixRtV_nxt  = 0.0d0
    xza_MixRtC_pre  = 0.0d0
    xza_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
    
    !---------------------------------------------------------------------
    ! $B<>=aK0OBD4@aK!$N<B9T(B
    !   $B%k!<%W$r2s$9$N$O(B, $B1@$K$D$$$F$@$1(B.  
    !---------------------------------------------------------------------
    LoopSvapPress: do s = 1, LoopNum
      
      !$B<>=aK0OBK!$G$O05NO$OJQ2=$5$;$J$$(B. 
      xz_ExnerAll = xz_Exner + xz_ExnerBasicZ
          
      !$B:#$^$G$KF@$i$l$?>qMp@.J,$NCM$r;CDjNL$H$_$J$9(B. $BE:$(;z$rDI2C(B
      xza_MixRtV_pre  = xza_MixRt(:,:,GasNum(s))   + xza_MixRtBasicZ(:,:,GasNum(s))
      xza_MixRtC_pre  = xza_MixRt(:,:,CloudNum(s)) + xza_MixRtBasicZ(:,:,CloudNum(s))
      xz_PotTemp_pre = xz_PotTemp
          
      Adjusting: do i = 0, ItrNum
        !---------------------------------------------------------------
        ! $BK0OB>x5$05$+$iK0OB:.9gHf$r5a$a$k(B
        !---------------------------------------------------------------
        !$B29EY(B
        xz_TempAll = ( xz_PotTemp_pre + xz_PotTempBasicZ ) * xz_ExnerAll

        !$BK0OB>x5$05$+$iK0OB:.9gHf$r7W;;(B($B4pK\>l$+$i$N:9(B). 
        xz_MixRtSat =                                                    &
          & xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll)              &
          &  * MolWtWet(CloudNum(s))                                     &
          &  / (MolWtDry * PressSfc * (xz_ExnerAll ** (CpDry / GasRDry)) )
        
        !$B5,3J2=$5$l$?@xG.(B
        xz_Gamma = xz_LatentHeatPerMass(SpcWetID(CloudNum(s)), xz_TempAll) / (xz_ExnerAll * CpDry)
        
        !$B6E7kNL$r5a$a$k(B. 
        !  $B6E7k$,@8$8$k>l9g$K$O(B, xza_MixRtV_pre - xz_MixRtSat $B$OI,$:@5$NCM$H$J$k(B.
        !  $B>xH/$,@8$8$k>l9g$K$O(B, $B>xH/NL$O(B - MixRtC $B$rD6$($k$3$H$O$J$$(B. 
        !  merge $B4X?t$rMQ$$$F>xH/NL$N>e8B$rDj$a$F$$$k(B. 
        xz_DelMixRt = &
          & merge( &
          &    (xza_MixRtV_pre - xz_MixRtSat), &
          &    - xza_MixRtC_pre, &
          &    (xza_MixRtV_pre - xz_MixRtSat) > - xza_MixRtC_pre &
          &   ) &
          &   / (1.0d0 + xz_Gamma * xz_DMixRtSatDPotTemp( &
          &           SpcWetID(CloudNum(s)), MolWtWet(CloudNum(s)), xz_TempAll, xz_ExnerAll & 
          &        ) )
        
        !$B$h$j??$K6a$$CM$r7W;;(B
        xz_PotTemp_nxt = xz_PotTemp_pre + xz_Gamma * xz_DelMixRt
        xza_MixRtV_nxt  = xza_MixRtV_pre  - xz_DelMixRt
        xza_MixRtC_nxt  = xza_MixRtC_pre  + xz_DelMixRt
        
        !$B7+$jJV$7$N$?$a$NJQ?tDj5A(B
        xz_PotTemp_pre = xz_PotTemp_nxt
        xza_MixRtV_pre  = xza_MixRtV_nxt
        xza_MixRtC_pre  = xza_MixRtC_nxt
        
      end do Adjusting
      
      !$BH?I|2~NI8e$NCM$r??$NCM$H$7$FJ]4I(B
      xz_PotTemp                = xz_PotTemp_nxt                 
      xza_MixRt(:,:,GasNum(s))   = xza_MixRtV_nxt - xza_MixRtBasicZ(:,:,GasNum(s))
      xza_MixRt(:,:,CloudNum(s)) = xza_MixRtC_nxt - xza_MixRtBasicZ(:,:,CloudNum(s))
      
    end do LoopSvapPress
    
  end subroutine MoistAdjustSvapPress2

  
!!!--------------------------------------------------------------------------!!!
  subroutine MoistAdjustNH4SH(xz_Exner, xz_PotTemp, xza_MixRt ) 
    !
    ! $B2=3XH?1~$N05J?9UDj?t(B Kp $B$rMQ$$$?K0OB<>=aD4@aK!(B
    ! $B$?$@$7(B, $B8=:_$O(B NH3 + H2S --> NH4SH $B$N@8@.H?1~$K$N$_BP1~(B
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none
    
    !$BJQ?tDj5A(B
    real(8),intent(in)   :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B%(%/%9%J!<4X?t(B
    real(8),intent(inout):: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B290L(B
    real(8),intent(inout):: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                          !$B6E=L@.J,$N:.9gHf(B
    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) :: PressAll
    real(8) :: TempAll
    real(8) :: Gamma
    real(8) :: DelMixRt
    
    integer :: i, k, num          ! $BE:$(;z(B  
            
    !---------------------------------------------------------------------
    ! $B=i4|2=(B
    !---------------------------------------------------------------------
    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
    PressAll  = 0.0d0
    TempAll   = 0.0d0
    Gamma     = 0.0d0
    DelMixRt  = 0.0d0
    num       = 0
  
    do k = RegZMin, RegZMax
      do i = RegXMin, RegXMax

        !$B<>=aK0OBK!$G$O05NO$OJQ2=$5$;$J$$(B. 
        ExnerAll = xz_Exner(i,k) + xz_ExnerBasicZ(i,k)
        PressAll = PressSfc * (ExnerAll ** (CpDry / GasRDry))        

        !$B:#$^$G$KF@$i$l$?>qMp@.J,$NCM$r;CDjNL$H$_$J$9(B. $BE:$(;z$rDI2C(B
        MixRtNH3_pre    = xza_MixRt(i,k,NH3Num)   + xza_MixRtBasicZ(i,k,NH3Num)
        MixRtH2S_pre    = xza_MixRt(i,k,H2SNum)   + xza_MixRtBasicZ(i,k,H2SNum)
        MixRtNH4SH_pre  = xza_MixRt(i,k,NH4SHNum) + xza_MixRtBasicZ(i,k,NH4SHNum)
        PotTemp_pre     = xz_PotTemp(i,k)
        
        !$B%k!<%W2s?t$N=i4|2=(B
        num = 0
        
        AdjustNH4SH: do
          !---------------------------------------------------------------
          ! $BJQ?t$N=i4|2=(B
          !---------------------------------------------------------------
          !$B29EY(B
          TempAll = (PotTemp_pre + xz_PotTempBasicZ(i,k)) * ExnerAll
          
          !$B5,3J2=$5$l$?H?1~G.(B (NH4SH $B#1(Bkg $B$KBP$9$kG.NL(B)
          Gamma = ReactHeatNH4SHPerMass / (ExnerAll * CpDry)
          
          DelMixRt = DelMixRtNH4SH(TempAll, PressAll,                         &
            &                      MixRtNH3_pre, MixRtH2S_pre, MixRtNH4SH_pre,             &
            &                      MolWtWet(NH3Num), MolWtWet(H2SNum), MolWtWet(NH4SHNum) )  

          !---------------------------------------------------------------
          ! $B$h$j??$K6a$$CM$r5a$a$kK0OB>x5$05$+$iK0OB:.9gHf$r5a$a$k(B
          !---------------------------------------------------------------
          ! NH4SH $B$N:.9gHf$r=$@5(B
          MixRtNH4SH_nxt  = MixRtNH4SH_pre + DelMixRt
          
          ! DelPress $B$r85$K(B, NH3 $B$H(B H2S $B$N:.9gHf$r=$@5(B
          MixRtNH3_nxt = MixRtNH3_pre - DelMixRt * MolWtWet(NH3Num) / MolWtWet(NH4SHNum)
          MixRtH2S_nxt = MixRtH2S_pre - DelMixRt * MolWtWet(H2SNum) / MolWtWet(NH4SHNum)
          
          !$B290L$r=$@5(B
          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, xza_MixRt(i,k,H2SNum), xza_MixRtBasicZ(i,k,H2SNum)
!            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
          
          !---------------------------------------------------------------
          ! $B<}B+H=Dj(B
          !---------------------------------------------------------------
          !$B290L$,<}B+$7$F$$$?$i7+$jJV$7$rCf;_$9$k(B
          if (abs((PotTemp_nxt - PotTemp_pre) / PotTemp_nxt) <= MinERROR   &
            & .OR. abs(PotTemp_nxt - PotTemp_pre) <= MinERROR             ) then 
            xz_PotTemp(i,k)         = PotTemp_nxt                 
            xza_MixRt(i,k,NH3Num)   = MixRtNH3_nxt - xza_MixRtBasicZ(i,k,NH3Num)
            xza_MixRt(i,k,H2SNum)   = MixRtH2S_nxt - xza_MixRtBasicZ(i,k,H2SNum)
            xza_MixRt(i,k,NH4SHNum) = MixRtNH4SH_nxt - xza_MixRtBasicZ(i,k,NH4SHNum)
            exit AdjustNH4SH
          end if
          
          !$B%k!<%W2s?t$,2s$j$9$.$?$i%(%i!<$rEG$$$FB3$1$k(B
          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                 
            xza_MixRt(i,k,NH3Num)   = MixRtNH3_nxt - xza_MixRtBasicZ(i,k,NH3Num)
            xza_MixRt(i,k,H2SNum)   = MixRtH2S_nxt - xza_MixRtBasicZ(i,k,H2SNum)
            xza_MixRt(i,k,NH4SHNum) = MixRtNH4SH_nxt - xza_MixRtBasicZ(i,k,NH4SHNum)
            exit AdjustNH4SH
          end if
          
          !$B7+$jJV$7$N$?$a$NJQ?tDj5A(B
          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, xza_MixRt ) 
    !
    ! NH3 + H2S --> NH4SH $B$N@8@.H?1~$N05J?9UDj?t(B Kp $B$rMQ$$$?K0OB<>=aD4@aK!(B
    !
    
    !$B0EL[$N7?@k8@6X;_(B
    implicit none
    
    !$BJQ?tDj5A(B
    real(8),intent(in)   :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B%(%/%9%J!<4X?t(B
    real(8),intent(inout):: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) 
                                                          !$B290L(B
    real(8),intent(inout):: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
                                                          !$B6E=L@.J,$N:.9gHf(B
    real(8) :: xz_PotTemp_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xz_PotTemp_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xza_MixRtNH3_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xza_MixRtNH3_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xza_MixRtH2S_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xza_MixRtH2S_nxt(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xza_MixRtNH4SH_pre(DimXMin:DimXMax, DimZMin:DimZMax) 
    real(8) :: xza_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) 
    
    integer            :: i
    integer, parameter :: ItrNum = 10
            
    !---------------------------------------------------------------------
    ! $B=i4|2=(B
    !---------------------------------------------------------------------
    xz_PotTemp_nxt = 0.0d0
    xza_MixRtNH3_nxt   = 0.0d0
    xza_MixRtH2S_nxt   = 0.0d0
    xza_MixRtNH4SH_nxt = 0.0d0
    
    xz_Gamma     = 0.0d0
    xz_DelMixRt  = 0.0d0

    !$B<>=aK0OBK!$G$O05NO$OJQ2=$5$;$J$$(B. 
    xz_ExnerAll = xz_Exner + xz_ExnerBasicZ
    xz_PressAll = PressSfc * (xz_ExnerAll ** (CpDry / GasRDry))
        
    !$B:#$^$G$KF@$i$l$?>qMp@.J,$NCM$r;CDjNL$H$_$J$9(B. $BE:$(;z$rDI2C(B
    xza_MixRtNH3_pre(:,:)    = xza_MixRt(:,:,NH3Num)   + xza_MixRtBasicZ(:,:,NH3Num)
    xza_MixRtH2S_pre(:,:)    = xza_MixRt(:,:,H2SNum)   + xza_MixRtBasicZ(:,:,H2SNum)
    xza_MixRtNH4SH_pre(:,:)  = xza_MixRt(:,:,NH4SHNum) + xza_MixRtBasicZ(:,:,NH4SHNum)
    xz_PotTemp_pre          = xz_PotTemp
    
    AdjustNH4SH: do i = 1, ItrNum
      !---------------------------------------------------------------
      ! $BJQ?t$N=i4|2=(B
      !---------------------------------------------------------------
      !$B29EY(B
      xz_TempAll = ( xz_PotTemp_pre + xz_PotTempBasicZ ) * xz_ExnerAll
      
      !$B5,3J2=$5$l$?H?1~G.(B (NH4SH $B#1(Bkg $B$KBP$9$kG.NL(B)
      xz_Gamma = ReactHeatNH4SHPerMass / ( xz_ExnerAll * CpDry )

      !NH4SH $B$N@8@.NL(B
      xz_DelMixRt = xz_DelMixRtNH4SH(xz_TempAll, xz_PressAll,           &
        &                      xza_MixRtNH3_pre, xza_MixRtH2S_pre, xza_MixRtNH4SH_pre,     &
        &                      MolWtWet(NH3Num), MolWtWet(H2SNum), MolWtWet(NH4SHNum) )  

      !---------------------------------------------------------------
      ! $B$h$j??$K6a$$CM$r5a$a$kK0OB>x5$05$+$iK0OB:.9gHf$r5a$a$k(B
      !---------------------------------------------------------------
      ! NH4SH $B$N:.9gHf$r=$@5(B
      xza_MixRtNH4SH_nxt  = xza_MixRtNH4SH_pre + xz_DelMixRt
      
      ! DelPress $B$r85$K(B, NH3 $B$H(B H2S $B$N:.9gHf$r=$@5(B
      xza_MixRtNH3_nxt = xza_MixRtNH3_pre - xz_DelMixRt * MolWtWet(NH3Num) / MolWtWet(NH4SHNum)
      xza_MixRtH2S_nxt = xza_MixRtH2S_pre - xz_DelMixRt * MolWtWet(H2SNum) / MolWtWet(NH4SHNum)
          
      !$B290L$r=$@5(B
      xz_PotTemp_nxt = xz_PotTemp_pre + xz_Gamma * xz_DelMixRt
      
      !$B%k!<%W$r2s$9$?$a$NJQ?tJQ2=(B
      xz_PotTemp_pre    = xz_PotTemp_nxt
      xza_MixRtNH3_pre   = xza_MixRtNH3_nxt 
      xza_MixRtH2S_pre   = xza_MixRtH2S_nxt 
      xza_MixRtNH4SH_pre = xza_MixRtNH4SH_nxt 

    end do AdjustNH4SH

    xz_PotTemp             = xz_PotTemp_nxt                 
    xza_MixRt(:,:,NH3Num)   = xza_MixRtNH3_nxt - xza_MixRtBasicZ(:,:,NH3Num)
    xza_MixRt(:,:,H2SNum)   = xza_MixRtH2S_nxt - xza_MixRtBasicZ(:,:,H2SNum)
    xza_MixRt(:,:,NH4SHNum) = xza_MixRtNH4SH_nxt - xza_MixRtBasicZ(:,:,NH4SHNum)
    
  end subroutine MoistAdjustNH4SH2
  
end Module MoistAdjust
